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

乘法逆元

在中国剩余定理的计算里,需要求一个数字在一个模下的逆元,也就是对于给定的 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 用引用的形式传入,返回值是最大公约数。


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

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

有句古话说得好:一时迭代一时爽,一直迭代一直爽。

这话我没说过——鲁迅.jpg
鲁迅:这锅我不背
这点困难怎么能阻挡我们寻找迭代算法的脚步!

注意到,递归的每一步都是关于 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,表示乘法逆元不存在。