python/TangDou/AcWing/NumberTheory/扩展欧几里得算法.md

扩展欧几里得算法

我们知道，计算两个正整数的最大公约数有两个比较常用的方法：更相减损术和辗转相除法，其中辗转相除法也叫欧几里得算法：gcd(a,b)=gcd(b,a \%b)

更相减损术和辗转相除法的主要区别在于前者所使用的运算是“减”，后者是“除”。从算法思想上看，两者并没有本质上的区别，但是在计算过程中，如果遇到一个数很大，另一个数比较小的情况，可能要进行很多次减法才能达到一次除法的效果，从而使得算法的时间复杂度退化为O(N)，其中N是原先的两个数中较大的一个。相比之下，辗转相除法的时间复杂度稳定于O(logN)。

一、欧几里得和扩展欧几里得

• 欧几里得算法，就是人们常说的辗转相除法，比较好理解，主要作用是求两个数最大公约数，最小公倍数也可方便的求出

• 扩展欧几里得就非常神奇了，主要作用是解不定方程， 即 a * x + b * y = c ，我们都知道可以有解，但不是唯一解

• 用 exgcd 可以很快求出 a * x + b * y = gcd(a,b) 的一个特解。如果 gcd(a,b) | c 即，有整数解，否则无整数解。

利用二元一次不定方程的通解形式一文中的证明，知道通解：

\large  x=x_0*k+\frac{b}{(a,b)}\cdot t 

\large  y=y_0*k-\frac{a}{(a,b)}\cdot t

其中t为任意整数,\large k=\frac{c}{gcd(a,b)}

二、扩展欧几里得实现过程

裴蜀等式求解过程:

应用欧几里得算法是一个递归的过程，将规模较大的问题转化为等价的小问题去解决，核心代码就是gcd(a,b)=gcd(b,a\%b),证明也略去了（笔者数学一般，见谅）。

它的边界情况出现在 gcd(k,0) 的时候，这时k就是 gcd(a,b) 啦。而扩展欧几里得算法则多了几个额外的步骤，在递归返回的时候做了一些小操作，使得我们可以得到 ax+by=gcd(a,b) 的一组解。思路如下：

1. 第一层： ax+by=gcd(a,b)，递归到下一层，gcd(a,b)=gcd(b,a\%b), 同时向下一层传 x 和 y ，设为x′和y′。
   
2. 第二层中：bx′+a\%b∗y′=gcd(b,a\%b) ,我们可以把 a\%b 表示为 a−⌊\frac{a}{b}⌋\cdot b ,代入整理可得

bx′+(a-⌊\frac{a}{b}⌋\cdot b)*y′=gcd(b,a\%b)

\Rightarrow 

ay′+b(x′−⌊\frac{a}{b}⌋ y′)=gcd(b,a\%b)=gcd(a,b)

比较可得，第一层中的 x 和 y 和第二层中 x′和 y′ 的关系为

\large x=y′,y=x′−⌊\frac{a}{b}⌋y′ 

3. 第二层的 $x′$和 $y′$ 又可以由第三层的 $x′′$和 $y′′$ 代入同样的求得，在递归终点的时，我们只需让 $x=1$ $y=0$即可。

注意因为上一层需要知道下一层的 x 和 y ，所以我们可以用传引用的方法，在递归返回时 x 和 y 就是下一层的 x′ 和 y′ 了。

扩展欧几里得解不定方程的步骤

#include <bits/stdc++.h>

using namespace std;
typedef long long LL;
//扩展欧几里得
LL exgcd(LL a, LL b, LL &x, LL &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    LL d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}
int main() {
    /* 原始方程 4x+6y=10 求x,y的一组解

    分析及求解步骤：
    此题中a=4,b=6,c=10
    1. 求d=gcd(a,b)=gcd(4,6)=2
    2. 判断 2|c即 10%2=0，方程有解，如果不能整除，此处可以直接返回方程无整数解
    3. 化简方程为 a/d*x+b/d*x=c/d 得到 4/2*x+6/2*y=10/2
       => 2x+3y=5 此方程与原始方程解是一样的，我们最终对这个方程下手。
    4. 如可解此方程呢？我们有一个好用的定理：扩展欧几里得可以解类似的方程
       ax+by=gcd(a,b)!
       但问题是按扩欧的描述，我们现在能计算的是 2x+3y=1的方程解！！！
       不是原方程2x+3y=5的解!!
       该怎么办呢？
       设 x0,y0为2x+3y=1的解，那么原方程的一组解就是：
       x1=5*x0,y1=5*y0!!!!

       为什么呢？原因很简单： 2x+3y=1 左右都乘以5
       2*(x*5)+3*(y*5)=5
       对比一下原方程 x1=x0*5,y1=y0*5啊~
    */
    LL x, y;
    exgcd(2, 3, x, y);
    printf("2x+3y=1 -> x=%lld y=%lld\n", x, y);
    printf("2x+3y=5 -> x=%lld y=%lld", 5 * x, 5 * y);
    return 0;
}

练习题I

P1082 [NOIP2012 提高组] 同余方程

#include <bits/stdc++.h>

using namespace std;
typedef long long LL;

// P1082 [NOIP2012 提高组] 同余方程

//扩展欧几里得
LL exgcd(LL a, LL b, LL &x, LL &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    LL d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

int main() {
    LL a, b;
    scanf("%lld %lld", &a, &b);
    LL x, y;
    exgcd(a, b, x, y);
    printf("%d\n", (x % b + b) % b);
    return 0;
}

练习题II

P5656 【模板】二元一次不定方程 (exgcd)

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;

//扩展欧几里得
LL exgcd(LL a, LL b, LL &x, LL &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    LL d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

int main() {
    int t;
    scanf("%d", &t);
    while (t--) {
        LL a, b, c;
        scanf("%lld %lld %lld", &a, &b, &c);
        LL x, y;
        LL d = __gcd(a, b); //使用库函数__gcd计算最大公约数d
        if (c % d)          //裴蜀定理,当gcd(a,b)不能整除c时，方程无解,输出-1
            puts("-1");
        else {
            //这步很值得总结：一个方程拿到手，第一步是化简，就是把a,b,c中都除以gcd(a,b),实现化简
            //举个栗子：2x+4y=8 -> x+2y=4 很直白吧，套路吧~ 此方程与原方程的角系一样吧~
            //当然，也可能化简不了，比如： 2x+3y=5,因为gcd(2,3)=1,化不化简都一样
            //但这里有一个事实，因为 a'=a/gcd(a,b),b'=b/gcd(a,b)后，a',b'肯定是互质。
            //令d'=gcd(a',b') 现在d'=1
            a /= d, b /= d, c /= d;

            //扩展欧几里得，一通用形式返回d,在有的题解中，不采用exgcd返回d,而是用__gcd计算d.
            //把扩展欧几里得算法定位为取特解x,y专用，也是比较直观简便的。
            //性能上嘛，由于辗转相除法不断的变更除数，速度上很快，两次也没有问题。
            //扩欧的形式：ax+by=gcd(a,b),求方程的一组解 x,y。现在的情况下，gcd(a,b)=1了，可以
            exgcd(a, b, x, y);
            x *= c, y *= c;         //结果(x,y)需要乘以c才能还原为ax+by=c的方程解。
            /*
            通解公式： x=x0+b/d*t
                      y=y0-a/d*t

            此时d=1,则通解公式：
                    x=x0+b*t
                    y=y0-a*t

            想要xmin => x0固定值，b固定值，只能变化的是t,每变化一次就增加、减少一个b
            */
            LL xmin = (x % b + b) % b; //利用最小正整数解的公式
            if (xmin == 0) xmin = b;   // xmin如果是0，需要特判,加一个b搞定~

            LL ymin = (y % a + a) % a;
            if (ymin == 0) ymin = a; // ymin如果是0，需要特判,加一个a搞定~

            // x最大值怎么求呢？简单：ax+by=c => ax=c-by => x = (c-by)/a
            // x要想最大，则y必须取最小值
            LL xmax = (c - b * ymin) / a;
            // y要想最大，则x必须取最小值
            LL ymax = (c - a * xmin) / b;

            if (xmax <= 0 || ymax <= 0)            //若最大值还不是正整数，则无正整数解
                printf("%lld %lld\n", xmin, ymin); //输出最小值
            else {
                //解的个数是最大 x 最小 x 的差除以模数 b 再加上1。
                LL sum = (xmax - xmin) / b + 1; //不足b时，也有1个解，+b就增加一个解
                printf("%lld %lld %lld %lld %lld\n", sum, xmin, ymin, xmax, ymax);
            }
        }
    }
    return 0;
}