这篇文章上次修改于 613 天前，可能其部分内容已经发生变化，如有疑问可询问作者。

乘法逆元

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

根据同余的定义，有 $b \mid (a a^* - 1)$ ，也就是存在整数 $k$ 使得 $bk=a a^*-1$ 。移一下项，就得到了 $a a^*-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\bmod b,b)=(b, a \bmod b)$ 。注意到 $b > (a \bmod 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 \bmod b)=(b, a-b\lfloor \frac{a}{b} \rfloor)$ ，我们将由此得出与 x，y 相关的递归关系。

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

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

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

这样，要求解 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{cases} x=y'\\[2ex] y=x'-\lfloor \frac{a}{b} \rfloor y' \end{cases}$

然后变成矩阵形式：

$\begin{pmatrix}x\\y\end{pmatrix}=\begin{pmatrix}0&1\\1&-d_1\end{pmatrix}\begin{pmatrix}x'\\y'\end{pmatrix}$

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

然后就可以一直乘下去：

$\begin{pmatrix}x\\y\end{pmatrix}=\begin{pmatrix}0&1\\1&-d_1\end{pmatrix}\begin{pmatrix}0&1\\1&-d_2\end{pmatrix}\begin{pmatrix}0&1\\1&-d_3\end{pmatrix}\cdots\begin{pmatrix}0&1\\1&-d_n\end{pmatrix}\begin{pmatrix}1\\0\end{pmatrix}$

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

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

$\begin{pmatrix}m_1&m_2\\n_1&n_2\end{pmatrix}\begin{pmatrix}0&1\\1&-d_1\end{pmatrix}=\begin{pmatrix}m_2&m_1-d_t m_2\\n_2&n_1-d_t n_2\end{pmatrix}$

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

再看一看初始值，就是单位阵 $\begin{pmatrix}m_1&m_2\\n_1&n_2\end{pmatrix}=\begin{pmatrix}1&0\\0&1\end{pmatrix}$ 。

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

现在梳理一下思路：

1. 第一步是给矩阵赋初始值。 $\begin{pmatrix}x&m\\y&n\end{pmatrix}=\begin{pmatrix}1&0\\0&1\end{pmatrix}$
2. 然后开始迭代。每一次迭代时，计算 $d = \lfloor \frac{a}{b} \rfloor$ 。
3. 更新矩阵，把矩阵右乘迭代矩阵的结果赋值给原先的矩阵。 $\begin{pmatrix}x&m\\y&n\end{pmatrix}\leftarrow\begin{pmatrix}x&m\\y&n\end{pmatrix}\begin{pmatrix}0&1\\1&-d\end{pmatrix}=\begin{pmatrix}m&x-d m\\n&y-d n\end{pmatrix}$
4. 进行一般欧几里得算法的迭代。 $\begin{pmatrix}a\\b\end{pmatrix}\leftarrow\begin{pmatrix}b\\a\bmod b\end{pmatrix}$
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 \pmod{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，表示乘法逆元不存在。