欧几里德算法及扩展欧几里德

来源：互联网发布：淘宝买家秀暴露狂女编辑：程序博客网时间：2024/05/17 01:47

欧几里德算法
欧几里德算法又称辗转相除法，用于计算两个整数a,b的最大公约数。其计算原理依赖于下面的定理：

定理：gcd(a,b) = gcd(b,a mod b)

证明：a可以表示成a = kb + r，则r = a mod b
假设d是a,b的一个公约数，则有
d|a, d|b，而r = a - kb，因此d|r
因此d是(b,a mod b)的公约数

假设d 是(b,a mod b)的公约数，则
d | b , d |r ，但是a = kb +r
因此d也是(a,b)的公约数

因此(a,b)和(b,a mod b)的公约数是一样的，其最大公约数也必然相等，得证

可以递归或者迭代写出程序：

int gcd(int n,int m)
{
int max,min;
int r;
n>m ? (max=n,min=m):(max=m,min=n);
while(min)
{
  r=max%min;
  max=min;
  min=r;
}
return max;
}

int gcd_(int n, int m)
{
int max,min;
n>m ? (max=n,min=m):(max=m,min=n);
if(min==0)
return max;
return gcd_(min,max%min);
}

扩展欧几里德算法不但能计算(a，b）的最大公约数，而且能计算a模b及b模a的乘法逆元，扩展欧几里德算法是用来在已知a, b求解一组p，q使得p * a+q * b = Gcd(a, b) (解一定存在，根据数论中的相关定理)。扩展欧几里德常用在求解模线性方程及方程组中。

把这个实现和Gcd的递归实现相比，发现多了下面的x,y赋值过程，这就是扩展欧几里德算法的精髓。
可以这样思考:
对于a' = b, b' = a % b 而言，我们求得 x, y使得 a'x + b'y = Gcd(a', b')
由于b' = a % b = a - a / b * b (注：这里的/是程序设计语言中的除法)
那么可以得到:
a'x + b'y = Gcd(a', b') ===>
bx + (a - a / b * b)y = Gcd(a', b') = Gcd(a, b) ===>
ay +b(x - a / b*y) = Gcd(a, b)
因此对于a和b而言，他们的相对应的p，q分别是 y和(x-a/b*y)

扩展欧几里德解一般二元一次方程a * x + b * y = c

从上面的过程可以看到，ax + by = gcd(a ,b)一定有解，再看一般形式ax + by=c，只有当且仅当c是gcd(a, b)的整数倍时才有解，否则无解。因为如果我们两边同除gcd(a, b)，会出现这样一种情况a'x + b'y=c / gcd(a, b)，如果c不是gcd(a,b)的整数倍，那么左边是整数，右边是分数，显然无解。

至此，到网上搜很多博客都会发现这句话“ax+by=c的求解可以先求出ax + by = gcd(a, b)，然后将x y扩大c / gcd(a,b)倍就可以了“，【而我自己无法证明这句话的正确性，反而觉得这句话是错的，在前面已经得出ax
+ by = gcd(a, b)的解是：

x = x0 + b / gcd(a, b) * k
y = y0 - a / gcd(a, b) * k （k为任意整数）

如果将x，y扩大c / gcd(a, b)倍后就成了

x = x0 * c / gcd(a, b) + b * c * k / (gcd(a, b) ^ 2)
y = y0 * c / gcd(a, b) - a * c * k / (gcd(a, b) ^ 2)（k为任意整数）

而我解出的通解是

x = x0 * c / gcd(a, b) + b / gcd(a, b) * k
y = y0 * c / gcd(a, b) - a / gcd(a, b) * k （k为任意整数）

上面两组通解明显不等价，下面的一组通解包含了上面的一组通解，通过代入原方程ax + by = c验证，下面的一组通解才是正确的。】【前面这段可忽略】所以，在求通解的过程中，只需要在特解x0 y0上扩大c / gcd(a, b)倍就可以了。


/*      Author: bcegkmqsw*/#include<stdio.h>#include<string.h>#include<stdlib.h>typedef __int64 LL;LL exgcd(LL a, LL b, LL &x, LL &y){    LL d, t;    if(b == 0)    {        x = 1, y = 0;        return a;    }    d = exgcd(b, a % b, x, y);    t = x, x = y, y = t - (a / b) * x;    return d;}int main(){    LL x, y, m, n, L, d, t, a, b, c;    while(scanf("%I64d%I64d%I64d%I64d%I64d", &x, &y, &m, &n, &L) != EOF)    {        a = m - n, b = L, c = y - x;        if(a < 0) a = -a, c = -c; // a > 0 => d > 0 => b > 0        d = exgcd(a, b, x, y);        if(m == n || c % d != 0)            printf("Impossible\n");        else        {            c /= d, t = c * x;            printf("%I64d\n", (t % b + b) % b);        }    }    return 0;}

不定方程ax+by=c 现在终于到了本文重点解二元一次不定方程。看起来扩展Euclid算法是不定方程的一种特殊情况实际上呢不定方程却是用Euclid算法解的。对于不定方程ax+by=c设gcd(a,b)=d如果ax+by=c有解则d|c这也是许多奥数题的切入点。所以一旦d不是c的约数那么 ax+by=c一定无解。当d|c时先求出ax’+by’=d=gcd(a,b)的x'和y'由于已经有ax’+by’=d要求ax’+by’=c将整个式子同乘c/d倍即可。则x=x'*c/dy=y'*c/d。由上一段可知
只要ax+by=c有一个解它就有无数个解。 Euclid算法还可以求解同余方程ax≡b(mod m)及其最小x。这其实和不定方程ax+my=b没有区别。不定方程和同余方程一般都有范围限制这其实也很容易解决就不说了扩展欧几里得解一般二元一次方程最小的正整数解不定方程ax+by=c。

0 0