算法学习笔记（37）——扩展欧几里得算法

扩展欧几里得算法

扩展欧几里得算法

欧几里得算法/辗转相除法（Euclidean algorithm）

辗转相除法，又称欧几里得算法（英语：Euclidean algorithm），是求最大公约数的算法。

两个整数的最大公约数是能够同时整除它们的最大的正整数。辗转相除法基于如下原理：两个整数的最大公约数等于其中较小的数和两数相除余数的最大公约数。

例如，252和105的最大公约数是21
$252=21\times 12;105=21\times 5$ ；
因 $252-105=21\times (12-5)=147$ ，所以147和105的最大公约数也是21。
在这过程中，较大的数缩小了，所以继续进行同样的计算可以不断缩小这两个数直至余为零。
注意：0与任何数的最大公约数即为该数。
这时，所剩下的还没有变成零的数就是两数的最大公约数。

裴蜀定理（Bézout 定理）

两数的最大公约数可以用两数的整数倍相加来表示，如 $21=5\times 105+(-2)\times 252$

对任意整数 $a, b$ ，存在一对整数 $x,y$ ，满足 $ax+by=gcd(a, b)$

证明：

在欧几里得算法的最后一步，即 $b = 0$ 时，显然有一对整数 $x=1, y=0$ ，使得 $a*1 + 0*0 = gcd(a,0)$ 。

若 $b>0$ ，则 $gcd(a,b) = gcd(b, a \bmod b)$ ;
假设存在一对整数 $x,y$ ，满足 $b*x + (a \bmod b) * y = gcd(b, a \bmod b)$ ;
由于 $b*x + (a \bmod b) * y = b*x + (a - b \lfloor a/b \rfloor) * y = a * y + b * (x - \lfloor a/b \rfloor) * y$ ;
则令 $x' = y, y' = (x - \lfloor a/b \rfloor) y$ ;
即得到 $ax' + by' = gcd(a,b)$ ;
对欧几里得算法的递归过程应用数学归纳法，可知裴蜀定理成立。

证毕

扩展欧几里得算法（Extended Euclidean algorithm）

扩展欧几里得算法（英语：Extended Euclidean algorithm）是欧几里得算法（又叫辗转相除法）的扩展。已知整数a、b，扩展欧几里得算法可以在求得a、b的最大公约数的同时，能找到整数x、y（其中一个很可能是负数），使它们满足贝祖等式

a x + b y = g c d (a, b)

$ax+by=gcd(a, b)$

假设 $d = gcd(a, b)$ ，当求出一组x与y的值之后，即可求出所有的x与y的值：

a (x - \frac{b}{d} \cdot k) + b (y + \frac{a}{d} \cdot k) = a x + b y - \frac{a b}{d} \cdot k + \frac{b a}{d} \cdot k = a x + b y = g c d (a, b)

$a(x - \frac{b}{d}\cdot k) + b(y + \frac{a}{d}\cdot k) = ax + by - \frac{ab}{d}\cdot k + \frac{ba}{d}\cdot k = ax + by = gcd(a, b)$

故所有的x与y的值可以用已经求出的一组解 $(x_0, y_0)$ 表示为：

{\begin{cases} x = x_{0} - \frac{b}{d} \cdot k \\ y = y_{0} + \frac{a}{d} \cdot k \end{cases}

$\begin{cases} x = x_0 - \frac{b}{d}\cdot k \\ y = y_0 + \frac{a}{d}\cdot k \end{cases}$

// 求a与b的最大公约数，并求出满足裴蜀定理 ax+by=gcd(a,b) 的一组x与y的值
int exgcd(int a, int b, int &x, int &y)
{
    /* 
     * 确定边界条件，如果b==0，则gcd(a,b)=a,则a*1+b*0=a一定满足裴蜀定理
     * 这种边界情况在代码中是作为递归的最底层出现的
     * 这意味着当这一层解出来的时候，我们需要将其返回到上一层的解，直到返回至最开始的一层。
     */
    if (b == 0) { x = 1, y = 0; return a; }
    
    /*
     * 如果 b>0, gcd(a, b) = gcd(b, a%b)
     * b*x + (a % b)y = b*x + (a - a / b * b)*y = gcd(a, b)
     * a*y + b(x - a / b * y) = gcd(a, b)
     * x' = y, y' = x - a/b*y
     */
    int d = exgcd(b, a % b, x, y);
    int z = x, y = x, x = z - a / b * y;
    return d;
}

精简版：

int exgcd(int a, int b, int &x, int &y)
{
    if (b == 0) { x = 1, y = 0; return a; }
    // 在计算 b 与 a%b 的最大公约数时交换 x 与 y 的值，方便后续计算
    int d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

题目链接：AcWing 877. 扩展欧几里得算法

#include <iostream>

using namespace std;

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

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

求解线性同余方程

给定整数 $a, b, m$ ，求一个整数 $x$ 满足 $a*x \equiv b \pmod m$ ，或者给出无解。因为未知数的指数为 $1$ ，所以我们称之为一次同余方程，也称线性同余方程。

$a*x \equiv b \pmod m$ 等价于 $a*x-b$ 是 $m$ 的倍数，不妨设为 $-y$ 倍。则 $a*x = m*(-y) + b$ ，于是，该方程可以改写为 $a*x+m*y = b$ 。

则该问题转化为，给定 $a, m, b$ ，求出一组x与y的值，使得 $a*x+m*y = b$ 成立。
而裴蜀定理则为，非定 $a, m, gcd(a,m)$ ，求出一组x与y的值，使得 $a*x+m*y=gcd(a, m)$ 成立

根据裴蜀定理及其证明过程，线性同余方程有解当且仅当 $gcd(a, m)|b$ 。

所以求解过程如下，记 $d = gcd(a, m)$ ：

求出 $a*x + m*y = d$ 的一组特解 $x_0, y_0$ （扩展欧几里得算法）
由于 $d|b$ ，则令 $x_0, y_0$ 同时乘以 $b/d$ ，就得到了 $a*x + m*y = b$ 的一组特解 $\frac{b}{d}x_0, \frac{b}{d}y_0$

事实上， $a*x+m*y = b$ 的通解可以表示为：

{\begin{cases} x = \frac{b}{d} x_{0} + k \frac{m}{d} \\ y = \frac{b}{d} y_{0} + k \frac{a}{d} \end{cases}

$\begin{cases} x = \frac{b}{d} x_0 + k \frac{m}{d} \\ y = \frac{b}{d} y_0 + k \frac{a}{d} \\ \end{cases}$

其中 $ k \in \mathbb{Z} $，$ d = gcd(a, m) $，$ x_0, y_0$ 是 $ax + my = d$ 的一组特解。

题目链接：AcWing 878. 线性同余方程

#include <iostream>

using namespace std;

typedef long long LL;

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

int main()
{
    int n;
    scanf("%d", &n);
    
    while (n -- ) {
        int a, b, m;
        scanf("%d%d%d", &a, &b, &m);
        int x, y;
        int d = exgcd(a, m, x, y);
        // 线性同余方程有解的充要条件：b是d的倍数
        if (b % d == 0) printf("%d\n", (LL)b / d * x % m);
        else puts("impossible");
    }
    
    return 0;
}