扩展欧几里得算法
乘法逆元
在中国剩余定理的计算里，需要求一个数字在一个模下的逆元，也就是对于给定的 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 的值就计算完了。
现在梳理一下思路：
- 第一步是给矩阵赋初始值。 \begin{pmatrix}x&m\\y&n\end{pmatrix}=\begin{pmatrix}1&0\\0&1\end{pmatrix}
- 然后开始迭代。每一次迭代时，计算 d = \lfloor \frac{a}{b} \rfloor 。
- 更新矩阵，把矩阵右乘迭代矩阵的结果赋值给原先的矩阵。 \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}
- 进行一般欧几里得算法的迭代。 \begin{pmatrix}a\\b\end{pmatrix}\leftarrow\begin{pmatrix}b\\a\bmod b\end{pmatrix}
- 当 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，表示乘法逆元不存在。
（封面pid：72665414）