欧几里德算法及其扩展

欧几里德算法其实就是辗转相除法，求最大公因数。

具体做法是：用较大数除以较小数，再用所得的余数（第一余数）去除除数，再用所得的余数（第二余数）去除第一余数，如此反复，直到最后余数是0为止。最后的除数就是这两个数的最大公约数。

记gcd（a，b）为整数a和b的最大公因数。

public static int gcd(int m,int n){
    return n == 0 ? m :gcd(n,m%n);
}

欧几里德算法的扩展——裴蜀（贝祖）等式

（最下方有相关知识点的代码实现模板，请自行选择是否越过文字解释）

对任何整数a,b和它们的最大公约数d，关于未知数x和y的线性方程（裴蜀等式）：

ax+by=m（m为d的倍数）当且仅当gcd（a,b）=d时有整数解

裴蜀等式有解时必然有无穷多个整数解，每组解x、y都称裴蜀数

特别地，方程ax+by=1有整数解当且仅当整数a,b互素。

x，y可以通过扩展欧几里德求解：

欧几里德算法停止的状态是：a=gcd ,b=0 ,那么，a'x+b'y=gcd 此时x=1,y为任意数。

因为，这时候只要a=gcd的系数是1，那么只要b的系数是0或者其他值（无所谓是多少，反正任意数乘以0都等于0但是a 的系数一定是1），这时，我们可以得到：a*1+b*0=gcd。

当然这是最终状态，但是我们可以以此反推到最初状态。

假设当前我们要处理的是求出a和b的最大公约数，并求出x和y使得ax+by=gcd，而我们已经求出了下一个状态：b和a%b的最大公约数，并求出了一组X1和Y1使得:b*X1+(a%b)*Y1=gcd,那么这两个相邻的状态之间是否存在一种关系呢？

令a%b=k，==> a = b*(a/b)+k 即 a%b = a - (a/b)*b

那么: gcd=b*X1 + (a-(a/b)*b)*Y1

              =b*X1 + a*Y1 - (a/b)*b*Y1

              =a*Y1 + b*(X1 - a/b*Y1)

对比 a*Y1 + b*(X1 - a/b*Y1) 和 ax+by=gcd ，求一组x和y使得ax+by=gcd

由此得到递推公式：

        x=Y1

        y=X1-a/b*Y1

但是需要注意的是，根据这个递推所求的x是ax+by=d的解，而扩展欧几里德算法就是在求a,b的最大公约数d=gcd(a,b)的同时，求出贝祖等式ax+by=m的一个解（X0,Y0）

所以实际上，X0=x*m*（d^（-1））

通项：（t是一个倍数）

        x=X0 + (b/gcd) * t        （所有的x对b同模）        求得的x1、x2、x3、x4之间是k倍的（b/d）

        y=Y0 - (a/gcd）* t         (所有的y对a同模）        求得的y1、y2、y3、y4之间是k倍的（a/d)

相关知识点代码实现模板如下：

/**
*扩展欧几里德算法
*x=y1
*y=x1-a/b*y1
*/
public class Main{
    static long x;
    static long y;
    public static long gcd(long m,long n){
        return n == 0 ? m : gcd(n,m%n);
    }
    //最小公倍数
    public staic long min(long a,long b){
        return a*b/gcd(a,b);
    }
    /**
    *扩展欧几里德
    *调用完成后x y是ax+by=gcd（a，b）的解
    */
    public static long ext_gcd(long a,long b){
        if(b==0){
            x=1;
            y=0;
            return a;
        }
        long res=ext_gcd(b,a%b);
        long x1=x;
        x=y;
        y=x1-a/b*y;
        return res;
    }
    /**
    *线性方程
    *ax+by=m 当m是gcd（a，b）倍数时有解
    *等价于ax=m mod b
    */
    public static long linearEquation(long a,long b,long m)throws Exception{
        long d=ext_gcd(a,b);
        //m不是gcd（a，b）的倍数，这个方程无解
        if（m%d！=0）throw new Exception（"无解"）；
        long n=m/d;
        x*=n;
        y*=n;
        return d;
    }
    /**
    *求线性同余方程组
    *x=a1(%m1)
    * =a2(%m2)
    * =a3(%m3)
    *x=a1+m1y1
    *x=a2+m2y2
    *==>m1y1-m2y2=a2-a1 这是一个线性方程，可解出y1 linearEquation（m1，-m2，a2-a1） 
    *代回x=a1+m1y1，得：x0=a1+m1y1 --> x=x0+k*min(m1,m2),得一个新方程：
    *x=x0(mod min(m1,m2))
    */
    public static long linearEquationGroup(long[] a,long[] m)throws Exception{
        int len=a.length;
        if(len == 0 && a[0] == 0){
            return m[0];
        }
        for(int i=1;i<len;i++){
            //这里往前看是两个方程
            long a2_a1=a[i]-a[i-1];
            long d=linearEquation(m[i-1],-m[i],a2_a1);
            //现在的x是y1，用y1求得一个特解
            long x0=a[i-1]+m[i-1]*x;
            long min=m[i-1]*m[i]/d;
            a[i]=(x0%min+min)%min;//x0变成正数
            m[i]=min;
        }
        //合并完成之后，只有一个方程：x=a[len](%m[len])
        //long d=linearEquation(1,m[len-1],a[len-1]);
        return a[len-1]%m[len-1];
    }
    /**
    *求逆元
    *ax = 1 (%mo),gcd(a,mo)=1
    *ax+mo*y=1
    */
    public static long inverseElement(long a,long mo) throws Exception{
        long d=linearEquantion(a,mo,1);//ax+mo*y=1
        x=(x%mo+mo)%mo;//保证x>0
        return d;
    }
}

 如果想要得到x大于0的第一个解：

b/=d;    x = (x0%b+b)%b