现在的位置: 首页 > 综合 > 正文

欧几里德&扩展以及求解线性方程学习总结–附上poj1061解题报告

2017年02月22日 ⁄ 综合 ⁄ 共 2468字 ⁄ 字号 评论关闭

欧几里德算法

欧几里德就是辗转相除法,调用这个gcd(a,b)这个函数求解a,b的最大公约数

公式:

gcd(a,b)=gcd(b,a%b);并且gcd(a,b)=gcd(b,a)=gcd(-a,b)=gcd(|a|,|b|)

代码:

int gcd(int a,int b)//递归
{
    if(b==0)
        return a;
    return 
        gcd(b,a%b);
}

int gcd(int a,int b)//递归简化
{
    return b ? gcd(b,a%b) : a;
}

int gcd(int a, int b)//迭代
{
     while(b != 0)
     {
       int r = b;
       b = a % b;
       a = r;
     }
     return a;
}

扩展欧几里德:

对于不完全为 0 的非负整数 a,b,gcd(a,b)表示 a,b 的最大公约数,必然存在整
数对 x,y ,使得 gcd(a,b)=ax+by,而且这样的解,是有很多个的;
解法
我们按欧几里得算法结束时来考虑
1. 最后b=0,那么gcd(a,b) = a = ax+by,那么此时的x=1,y=0

2. b!=0的时候我们就按照欧几里得公式 gcd(a,b) = gcd(b,a%b)推导一般情况下x,y的值

gcd(a,b) = ax+by = gcd(b,a%b) = bx1+(a%b)y1;

那么a%b = a-a/b*b  注意这a/b 是语言除法,取整

那么就有bx1+(a%b)y1 = bx1+(a-a/b*b)y1 将其展开整理得 ax+by = ay1 + b(x1-a/b*y1) 

那么得到一般情况下的 x = y1, y = x1-a/b*y1  既x 和 y 的解是通过递归得到的一组解,递归到最后就是 x = 1,y = 0 那么实际的 gcd(a,b) = ax+by 的解就可以在递归中得到

代码如下

int exgcd(int a,int b,int &x,int &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return a;
    }
    int r=exgcd(b,a%b,x,y);
    int t=x;
    x=y;
    y=t-a/b*y;
    return r;
} //递归方法

int exgcd(int m,int n,int &x,int &y)
{
    int x1,y1,x0,y0;
    x0=1; y0=0;
    x1=0; y1=1;
    x=0; y=1;
    int r=m%n;
    int q=(m-r)/n;
    while(r)
    {
        x=x0-q*x1; y=y0-q*y1;
        x0=x1; y0=y1;
        x1=x; y1=y;
        m=n; n=r; r=m%n;
        q=(m-r)/n;
    }
    return n;
}//非递归方法

求解线性方程:
如果我们有线性方程: ax+by=c,(a,b,c一直)如果此方程有解,那么一定有 c%gcd(a,b) = 0,为什么呢?首先

a%gcd(a,b) = 0 = b%gcd(a,b); 那么c%gcd(a,b) = (ax+by)/gcd(a,b) = ax%gcd(a,b)+by%gcd(a,b) = 0+0=0

那么,接下来求解(x,y) ,首先我们通过exgcd(a,b) 就可以得到一组 x0,y0 ,那么我们现在通过这一组解得到通解

我们取另一组未知解(x,y);那么就有 <1式>ax0+by0=ax+by=gcd(a,b);//(不管哪一组解,他们的gcd,也就是最大公约数都是一样的)

那么我们根据<1式>
<2式>a(x-x0) = b(y0-y),现在我们定义 d = gcd(a,b);

<2式>我们同除 d 得
<3式>
a/d(x-x0) = b/d(y0-y);那么我们可以知道 a/d 和 b/d 现在互素那么根据<3式>两个互素的数分别乘一个数要相等,那么就有x-x0 = t(b/d)   y0-y = t(a/d);t为整数,也就是x-x0一定要是b/d的倍数,那么就可以得到 在有一组解 x0,y0 的情况下其他的解的通解如下式:

x=x0+t(b/d)

y=y0-t(a/d)

那么这样的解有无数多个,一般实际会要求一个条件下的最大或最小,

先给出求x最小的一个方法,下面引入poj1061青蛙的约会,就是要求所有的解当中的x最小

我们通过 exgcd(a,b) 求得了(x0,y0), 但是注意它不是满足ax+by=c 的一个解 ,因为我们只求得ax0+by0=gcd(a,b),但是这里c 和gcd(a,b) 还有一个倍数的关系,通过我们求得的ax0+by0=gcd(a,b) 两边同时 乘以 c/gcd(a,b) 才得到

ax0*(c/gcd(a,b)) + by0*(c/gcd(a,b)) = c;

那么根据上面的通解 x = x0+t(b/d);要使得x最小,那么就是取 t 最小 ,有一个技巧k=b/d,min = (x%k+k)%k;

poj1061

分析:那么假设跳了T次以后,青蛙A的坐标便是(x+m*T)%L,青蛙B的坐标为(y+n*T)%L,它们能够相遇的情况为(x+m*T)%L = (y+n*T)%L,那么他们之间跳的相对距离列式:(x+m*T) - (y+n*T)==P*L,其中P为某一个整数,变形一下得到(n-m)*T-P*L==x-y
我们设a=(n-m),b=L,c=x-y,T=x,P=y.于是便得到ax+by==c,直接利用欧几里得扩展定理可以得到一组x0,y0,接着用上面的方法即可:

上马:

// 132K 32MS
#include <stdio.h>

#define ll __int64

ll x,y,a,b,c,d;

ll exgcd(ll a,ll b)
{
    if(!b) {x = 1; y = 0; return a;}
    ll d = exgcd(b,a%b);
    ll X = x;
    x = y;
    y = X - a/b*y;
    return d;
}

int main()
{
    ll X,Y,n,m,L;
    while(~scanf("%I64d%I64d%I64d%I64d%I64d",&X,&Y,&m,&n,&L))
    {
        a = n-m; b = L; c = X-Y;
        d = exgcd(a,b);//先得到最大公约数,同时得到一组解x,y

        if(c%d != 0)//无解
        {
            printf("Impossible\n");continue;
        }

        ll k = b/d;//最小解
        x *= (c/d);//
        x = (x%k+k)%k;
        printf("%I64d\n",x);
    }
    return 0;
}

个人愚昧观点,欢迎指正和讨论

抱歉!评论已关闭.