扩展欧几里得算法

---

乘法逆元

在中国剩余定理的计算里，需要求一个数字在一个模下的逆元，也就是对于给定的 a，b，找到方程 $a{a}^{\ast }\equiv 1(\mathrm{mod}b)$ 的一个整数解 $a^{\ast }$ 。接下来我们分析一下这个方程背后隐藏着什么。

根据同余的定义，有 $b\mid (a{a}^{\ast }-1)$ ，也就是存在整数 $k$ 使得 $bk=a{a}^{\ast }-1$ 。移一下项，就得到了 $a{a}^{\ast }-bk=1$ 。

这个形式恰好符合裴蜀定理 $ax+by=1$ 的形式，于是 $(a,b)=1$ ，这表明a，b互质是逆元存在的必要条件。同样可以证明：a，b互质是 a 在模 b 下存在逆元的充分条件。

---

欧几里得算法

从上面的论述可以看出，要求一个数的逆元，实际上就是找到裴蜀定理 $ax+by=(a,b)=1$ 中的 x 和 y。实际上这一步工作可以通过辗转相除法完成。

碾转相除法的依据是一个基本的事实： $(a,b)=(a-kb,b),k\in \mathbb{Z}$ 。

这样，在 a > b时，令 $k=\lfloor \frac{a}{b}\rfloor$ ，我们得到 $(a,b)=(a-b\lfloor \frac{a}{b}\rfloor ,b)=(a\mathrm{mod}b,b)=(b,a\mathrm{mod}b)$ 。注意到 $b>(a\mathrm{mod}b)$ ，所以这个操作可以继续进行下去，直到 b=0，此时 $(a,0)=a$ 。

这种求最大公约数的算法叫做辗转相除法，又叫欧几里得算法。

根据上面的描述，很容易得出一个碾转相除法的递归实现：

int gcd(int a, int b) {
    if(a < b) return gcd(b, a);
    if(b == 0) return a;
    return gcd(b, a % b);
}
 

对这个函数进行尾递归内联优化，就得到了它的迭代版本：

int gcd(int a, int b) {
    if(a < b) return gcd(b, a);
    while(b != 0) {
        int t = a % b;
        a = b; b = t;
    }
    return a;
}
 

辗转相除法是计算最大公约数很快的一种方式，但是它的潜力远不止此。下面我们要用它完成更加神奇的事情。

---

扩展欧几里得算法

欧几里得算法的另一个用途就是求解 $(a,b)=ax+by$ 中的 x 和 y。

在欧几里得算法里，递归部分的核心是这样的： $(a,b)=(b,a\mathrm{mod}b)=(b,a-b\lfloor \frac{a}{b}\rfloor )$ ，我们将由此得出与 x，y 相关的递归关系。

根据裴蜀定理，我们可以找到四个整数 $x,y,x^{\prime },y^{\prime }$ ，使得 $ax+by=(a,b)$ 且 $b{x}^{\prime }+(a-b\lfloor \frac{a}{b}\rfloor )y^{\prime }=(b,a-b\lfloor \frac{a}{b}\rfloor )$ 。

这两个等式的右侧是相等的，于是我们得到：$ax+by=b{x}^{\prime }+(a-b\lfloor \frac{a}{b}\rfloor )y^{\prime }$ ，整理一下就是 $a(x-y^{\prime })+b(y-(x^{\prime }-\lfloor \frac{a}{b}\rfloor y^{\prime }))=0$ 。

我们希望这个等式对一切 a，b 都成立，于是 $x=y^{\prime }$ 且 $y=x^{\prime }-\lfloor \frac{a}{b}\rfloor y^{\prime }$ 。

这样，要求解 x 和 y，只需要求解 x’，y’，由于 x’，y’ 对应的问题规模更小，所以可以进行递归的运算。最后找一下边界条件：当 b=0 时， $(a,0)$ 对应 $x=1$ ， $y=0$ 。

这种算法的思想和欧几里得算法是一致的，而且它在求出最大公约数的同时，求解了一组适于裴蜀定理的系数，所以叫做扩展欧几里得算法。

递归算法也是很好写的：

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

x 和 y 用引用的形式传入，返回值是最大公约数。

---

扩展欧几里得算法的迭代形式

有一个不太令人高兴的事实：扩展欧几里得算法的递归是在函数中央的，不能进行尾递归内联。

有句古话说得好：一时迭代一时爽，一直迭代一直爽。

鲁迅：这锅我不背 这点困难怎么能阻挡我们寻找迭代算法的脚步！

注意到，递归的每一步都是关于 x’，y’ 的线性组合，所以我们可以考虑用类似于线性递推数列的方式，用矩阵的形式改写递归式，从而利用矩阵乘法的结合律找到迭代算法。

这是我们找到的递归式：

$\{\begin{array}{l}x=y^{\prime }\\ y=x^{\prime }-\lfloor \frac{a}{b}\rfloor y^{\prime }\end{array}$

然后变成矩阵形式：

$\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}0& 1\\ 1& -d_1\end{array}\right)\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)$

这里 $d_1$ 代表第一次迭代时的 $\lfloor \frac{a}{b}\rfloor$ 。

然后就可以一直乘下去：

$\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}0& 1\\ 1& -d_1\end{array}\right)\left(\begin{array}{cc}0& 1\\ 1& -d_2\end{array}\right)\left(\begin{array}{cc}0& 1\\ 1& -d_3\end{array}\right)\cdots \left(\begin{array}{cc}0& 1\\ 1& -d_n\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)$

在每次迭代时，都可以从左向右计算一个矩阵，这样就达到了迭代的目的。

假设在一次迭代之前计算的矩阵之积为 $\left(\begin{array}{cc}{m}_1& m_2\\ n_1& n_2\end{array}\right)$ ，然后下一个矩阵 $\left(\begin{array}{cc}0& 1\\ 1& -d_t\end{array}\right)$ 是右乘上去的：

$\left(\begin{array}{cc}{m}_1& m_2\\ n_1& n_2\end{array}\right)\left(\begin{array}{cc}0& 1\\ 1& -d_1\end{array}\right)=\left(\begin{array}{cc}{m}_2& m_1-d_t{m}_2\\ n_2& n_1-d_t{n}_2\end{array}\right)$

把所有矩阵的积计算完成之后，得到一个最终的矩阵 $\left(\begin{array}{cc}{M}_1& M_2\\ N_1& N_2\end{array}\right)$ ，这个矩阵要满足 $\left(\begin{array}{c}x\\ y\end{array}\right)\left(\begin{array}{cc}{M}_1& M_2\\ N_1& N_2\end{array}\right)=\left(\begin{array}{c}1\\ 0\end{array}\right)$ ，展开整理，就得到 $x=M_1$ ， $y=N_1$ 。

再看一看初始值，就是单位阵 $\left(\begin{array}{cc}{m}_1& m_2\\ n_1& n_2\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$ 。

所以这里可以有一个小(负)优化，在迭代时直接把矩阵写成 $\left(\begin{array}{cc}x& m_2\\ y& n_2\end{array}\right)$ ，这样迭代完成之后 x 和 y 的值就计算完了。

现在梳理一下思路：

1. 第一步是给矩阵赋初始值。 $\left(\begin{array}{cc}x& m\\ y& n\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$
2. 然后开始迭代。每一次迭代时，计算 $d=\lfloor \frac{a}{b}\rfloor$ 。
3. 更新矩阵，把矩阵右乘迭代矩阵的结果赋值给原先的矩阵。 $\left(\begin{array}{cc}x& m\\ y& n\end{array}\right)\leftarrow \left(\begin{array}{cc}x& m\\ y& n\end{array}\right)\left(\begin{array}{cc}0& 1\\ 1& -d\end{array}\right)=\left(\begin{array}{cc}m& x-dm\\ n& y-dn\end{array}\right)$
4. 进行一般欧几里得算法的迭代。 $\left(\begin{array}{c}a\\ b\end{array}\right)\leftarrow \left(\begin{array}{c}b\\ a\mathrm{mod}b\end{array}\right)$
5. 当 b = 0 时，计算结束，此时 x 和 y 已经计算完成，并且最大公约数储存在 a 中。 给一个迭代版本的代码：

int exgcd(int a, int b, int& x, int& y) {
    if(a < b) return exgcd(b, a, y, x);
    int m = 0, n = 1;
    x = 1; y = 0;
    while(b != 0) {
        int d = a / b, t;
        t = m; m = x - d * t; x = t;
        t = n; n = y - d * t; y = t;
        t = a % b; a = b; b = t;
    }
    return a;
}
 

---

扩展欧几里得算法与乘法逆元

当我们用扩展欧几里得算法找到一组 x 和 y 满足 $ax+by=1$ 时，就可以得出 $ax\equiv 1(\mathrm{mod}b)$ ，也就是说 x 是 a 的乘法逆元。

显然，如果 x 是 a 的乘法逆元，那么所有的 $x+kb,k\in \mathbb{N}$ 也是 a 的乘法逆元。这表明，一定有 a 的一个乘法逆元在区间 $[0,b)$ 内。

这是找到 n 在模 p 下的乘法逆元的一个算法，使用了上面的 exgcd：

int inv(int n, int p) {
    int x, y;
    if(exgcd(n, p, x, y) == 1) {
        x = x % p;
        return x >= 0 ? x : p + x;
    } else {
        return -1;
    }
}
 

如果返回 -1，表示乘法逆元不存在。