求解最大公约数——欧几里得算法及其（解同余方程）拓展欧几里得

最大公约数的求法中最过著名的莫过于欧几里得辗展相除法，它有两种形式（递归与非递归，其实是一样的，任何递归都可以写成非递归），下面看看GCD的实现代码：

/***求a,b最大公约数***/
long long gcd(long long a, long long b) {
        if(b == 0)
                return a;
        else
                return gcd(b, a % b);
}

gcd(a,b) = gcd(b,a mod b) (不妨设a>b 且r=a mod b ,r不为0)

证法一

a可以表示成a = kb + r（a，b，k，r皆为正整数），则r = a mod b

假设d是a,b的一个公约数，记作d|a,d|b，即a和b都可以被d整除。

而r = a - kb，两边同时除以d，r/d=a/d-kb/d=m，等式左边可知m为整数，因此d|r

因此d也是（b,a mod b）的公约数

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

证法二

第一步：令c=gcd(a,b），则设a=mc，b=nc

第二步：可知r =a-kb=mc-knc=(m-kn)c

第三步：根据第二步结果可知c也是r的因数

第四步：可以断定m-kn与n互素【否则，可设m-kn=xd，n=yd，（d>1），

则m=kn+xd=kyd+xd=(ky+x)d，则a=mc=(ky+x)dc，b=nc=ycd，故a与b最大公约数≥cd，而非c，与前面结论矛盾】

从而可知gcd(b,r)=c，继而gcd(a,b)=gcd(b,r），得证

以上两种方法实质一样的。

---

扩展的欧几里德算法是求如a * x + b * y = (a, b) 这样的整数解的，可以仿照欧几里德算法得出答案。程序如下：

/***扩展的欧几里德算法a*x + b*y = Gcd(a,b)的一组整数解，结果存在x,y中***/
void extend_gcd(long long a, long long b, long long& x, long long &y) {
        if(b == 0) {
                x = 1;
                y = 0;
                return;
        }
        extend_gcd(b, a % b, x, y);
        long long tmp = x;
        x = y;
        y = tmp - a / b * y;
}

上述程序只是得到了一组解，很显然解是不唯一的：x增加b, y减少a一定是原方程的一组解：

a *  (x + b)  + b * (y - a) = a * x + b * y = (a, b)。

然而在应用上，往往并不是如此简单，很多时候会求解不定方程a * x + b * y = n。这个时候还是应用上面的算法：

1. 求(a，b), 设c = (a，b)，如果! c|n，则不存在整数解。因为将上式左右两边都除以c，可以知道，左边为整数，右边为非整数，故矛盾。
2. 将左右两边同时除以c,设得到新的方程为a' * x + b' * y = n'，应用上述算法求a' * x + b' * y = 1的解（由第一步知道(a'，b') = 1）。设结果为x', y'。
3. x = x' * n' , y = y' * n'是方程a * x + b * y = n。这个比较好理解，将a' * x + b' * y = 1两边同时扩大n'倍就行了。
4. x = x' * n' + t * b, y = y' * n' - t * a（t为整数)是原方程a * x + b * y = n的所有解。

线性同余方程实现代码：

/*  扩展欧几里得定理  扩展欧几里得定理：对于两个不全为0的整数a、b，必存在一组解x,y，  使得ax+by==gcd(a,b);  拓展欧几里得实现  下面面的程序中,x和y用全局变量保存  int gcd(int a,int b)  {    int t,d;    if(b==0)    {        x=1;        y=0; //此时b==0，也就是说gcd(a,0)==a。原式变为ax+by==a=gcd(a,b)--> x==1,y==0        return a; //返回a,b最大公约数的值     }    d=gcd(b,a%b); //欧几里得求最大公约数应用     t=x;    x=y;    y=t-(a/b)*y;  //不明处2    return d;     //返回a,b最大公约数的值  }   //不明处2 解释 ax+by==gcd(a,b)（1）     说明规则，x,y表示第一次递归时的值，x1,y1表示第二次递归时的值。那么gcd(a,b)==gcd(b,a%b)，同时都代入式1，有ax+by==b*x1+(a%b)*y1。将右边变形一下    b*x1+(a%b)*y1==b*x1+(a-(a/b)*b)*y1==a*y1+b*(x1-(a/b)*y1)，最终得到ax+by==a*y1+b*(x1-(a/b)*y1)    也就是说： 上一深度的x等于下一深度的y1，上一深度的y等于下一深度的x1-(a/b)*y1。    *需要注意，上面推导时用的除法都是整型除法    因此，得到了不定式ax+by==gcd(a,b)的一组解，x、y。    那么对于一般的不定式ax+by==c，它的解应该是什么呢。很简单，x1=x*(c/gcd(a,b)),y1=y*(c/gcd(a,b)) 再深入一点，就解出这么一组解其实一般来说是解决不了什么问题的，我们现在要得到所有的解，那么这所有的解究竟是什么呢？假设d=gcd(a,b). 那么x=x0+b/d*t; y=y0-a/d*t;其中t为任意常整数。  */ #include<stdio.h>#include<stdlib.h>/*  求关于 x 的同余方程 ax ≡ 1 (mod b)的最小正整数解。  即求ax=mb+r 1=nb+r  变形ax+(n-m)b=1，此方程即拓展欧几里德的应用ax+by=gcd(a,b),（n-m相当于y）  事实上ax ≡1(mod b) 有解的必要条件是gcd(a,b)|1，即gcd(a,b)=1；   使用拓展欧几里得可知 ax+by=1（x,y是整数） */int  Extragcd(int a,int b,int *x,int *y){int d,t;if(b==0)  //递归调用终止条件，当根据欧几里得辗转相除法则，余数为0停止 {*x=1;*y=0;return a;}else{d=Extragcd(b,a%b,x,y);t=*x;     //根据下一个x1,y1的值，倒推前一个x,y的值 *x=*y;*y=t-a/b*(*y);return d;}}int main(){int a,b,x,y;int ans;freopen("mod.in","r",stdin);freopen("mod.out","w",stdout); scanf("%d%d",&a,&b);ans=Extragcd(a,b,&x,&y);if(ans!=1) return 0;//根据若x是方程的一个解，则方程的所有解为x+k*b k为整数     x=x%b; //保证最小的正整数解x ,且x属于{0,1,2,3...b-1}     while(x<=0)        x+=b ;    printf("%d\n",x);    return 0;} 

ax+by=c问题

问题：ax+by=c，已知a、b、c，求解使该等式成立的一组x，y。其中a、b、c、x、y均为整数

a,b的最大公约数为gcd(a,b)。如果c不是gcd(a,b)的倍数，则该等式无解，因为等式左边除以gcd(a,b)是整数，

而等式右边除以gcd(a,b)后为小数。(根据解方程的时候，在等式的左右两边同时除以非0的整数，等式依然成立)

因此，只有当c是gcd(a,b)的倍数的时候，该等式有解。这样，可以通过计算使ax1+by1=gcd(a,b)成立的x1、y1，

然后有x=(c/gcd(a,b))*x1，y=(c/gcd(a,b))*y1，得到x,y。

问题就被转换为求使ax+by=gcd(a,b)成立的一组x,y。这可以用扩展欧几里德算法求解。如下：

根据欧几里得的辗转相除法，最终如果b为零，则gcd(a,b)=a，那么x=1,y=0为一组解。(一组重要的特解)

如果b不为零，根据欧几里德定理，可以设

ax1+by1=gcd(a,b)=gcd(b,a%b)=bx2+(a%b)y2=bx2+(a-(a/b)*b)y2

进一步化简可以得到 -> a*y2+b*(x2-(a/b)*y2)

化简后有x1=y2，y1=x2-(a/b)y2。因此x1,y1依赖于x2,y2，同理依次类推递归调用求出x3,y3,x4,y4……，

类似于辗转相除，直到b=0时，求出xn,yn，便可以由后往前倒推出x1,y1的值。

扩展欧几里德算法：

// 扩展欧几里德算法,解gcd(a, b) = ax + by// 结果存储在x,y中,用户调用时保证a、b、c都是整数// 返回a,b的最大公约数int extended_euclid(int a, int b, int &x, int &y){    if(b == 0) // gcd(a, b) == a    {        x = 1;        y = 0;        return a;    }    int n = extended_euclid(b, a%b, x, y);    int tmp = x;    x = y;    y = tmp - static_cast<int>(a/b)*y;    return n;}

用户调用linear_equation求解线性方程：

// 等式ax+by=c,已知a、b、c,求x和y。// 解该线性方程等同于解同余式ax = c(mod b)// 返回值表示是否有解,true有解,false无解bool linear_equation(int a, int b, int c, int &x, int &y){    int n = extended_euclid(a, b, x, y);    if(c%n)        return false;    int k = c/n;    x *= k;    y *= k;    return true;}

linear_equation函数也可以用来解同余式ax=c(mod b)。

由ax=c(mod b)，可以得到ax = mb+r；c = nb+r。化简可以得到ax+(n-m)b=c。调用linear_equation可以求出x。