乘法逆元

发表于 更新于 分类于

乘法逆元的定义及用途，以及求逆元的三种算法：扩展欧几里得、费马小定理、递推求逆元。

乘法逆元

首先，数学上的乘法逆元就是指直观的倒数，即 \(a\) 的逆元是 \(\frac{1}{a}\)，也即与 \(a\) 相乘得 1 的数。\(ax=1\)，则\(x\)是\(a\)的乘法逆元。

这里我们讨论关于取模运算的乘法逆元，即对于整数 \(a\)，与 \(a\) 互质的数 \(b\) 作为模数，当整数 \(x\) 满足 \(ax \bmod b \equiv 1\) 时，称 \(x\) 为 \(a\) 关于模 \(b\) 的逆元，代码表示就是a * x % b == 1。

在算法竞赛中，经常会遇到求解数据很大，则输出模 \(10^{9}+7\) 的解这类要求。加法、减法、乘法等操作，基于同余理论直接取模即可。但遇到除法时，某步中间结果不一定能完成整除，就无法求解了。

举个例子：求3 * 6 / 3 对 7 取模的结果。我们直接算出3 * 6 / 3的结果是6，对7取模得最终答案是 6 。

但我们通常面对的问题是中间结果超过int甚至long long 的范围，而不得不在每一步基于同余理论取模，我们用这个例子尝试一下：

还是求 3 * 6 / 3 % 7

第一步：3 * 6 == 18，18 % 7 == 4

第二步：4 这个中间结果再做 4 / 3 无法整除，就无法进行下去了。

但我们可以求出除数 3 关于模数7的逆元 5（根据逆元定义，5 符合 3 * 5 % 7 == 1），从而，用乘以5代替除以3。

上述第二步除法变乘法： 4 * 5 == 20，20 % 7 == 6

从而也计算出了正确的结果 6 。

乘法逆元的作用就是：

设 \(m\) 是一个很大的数，\(a\)、\(b\)已知，预期要计算（假设答案为 \(c\)）：

\[ m / a \bmod b \]

对于 \(a\) 的逆元 \(d\)，能够满足

\[ m \cdot d \bmod b = m / a \bmod b = c \]

在有些问题中，无法计算最终值很大的 \(m\) ，只能得到基于同余的一个中间值 \(m \bmod b = e\) 来计算 \(e / a \bmod b\) ，而 \(e\) 可能无法整除 \(a\)，就可以用 \(a\) 的逆元 \(d\)，来计算 \(e \cdot d \bmod b\)。

故而我们需要一个算法求 除数 的 取模逆元 ，从而在四则运算取模的任务中，用逆元将除法转为乘法。

方法1：扩展欧几里得

先暂时将逆元的事放一放，来看下扩展欧几里得。

首先大多入门选手知道求两个数最大公约数的算法，即辗转相除：

1
int GCD(int a, int b) {return b ? GCD(b, a % b) : a}

扩展欧几里得算法则是求 \(ax + by = GCD(a, b)\) 的一组可行解:

设两个式子

\[ ax + by=GCD(a, b) \]

\[ bx' + (a \bmod b)y'=GCD(b, a \bmod b) \]

由欧几里得算法知 \(GCD(a, b) = GCD(b, a \bmod b)\)

所以

\[ ax + by = bx' + (a \bmod b)y' \]

而

\[ a \bmod b = a - kb \]

其中 \(k = \left \lfloor a / b \right \rfloor\)，\(\lfloor \rfloor\)表示向下取整。

所以有

\[ ax + by = bx' + (a - kb)y' \]

展开移项得：

\[ ax + by = ay' + b(x'-ky') \]

根据系数对应关系，可以设

\(x = y'\)、\(y = x'-ky'\) 来进一步求一个可行解。

递归地用 \(x'\)、\(y'\) 表示“上一步”的\(x\)、\(y\)，就能递归地把问题转换成

\[ bx' + (a \bmod b)y'=GCD(b, a \bmod b) \]

类似 GCD() 的递归终点，当扩展欧几里得算法ExGCD()的\(a'x'+b'y'=GCD(a',b')\)中的\(b'\) 为 0 时，可以得到递归终点的 \(x'=1, y'=0\)，层层回溯套用前面等式

\[ x = y' \]

\[ y = x' - ky' \]

就能得到 \(ax+by=GCD(a, b)\) 的一组可行解。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
typedef long long LL;
LL ExGCD(LL a, LL b, LL &x, LL &y)
{
    // x, y 为引用传参，故最终程序结束后，x,y会被赋值为可行解
    if(b == 0)
    {
        // 递归终点，ax+by=GCD(a,b)的b为0，故方程变为
        // ax=a，则可行解可以是 x=1, y=0
        x = 1, y = 0;
        return a;
    }
    LL d = ExGCD(b, a % b, x, y), t = x;
    x = y, y = t - a / b * x;
    return d;  // 这里返回值是GCD(a,b)的结果，即最大公约数
}

扩展欧几里得求逆元

了解了扩展欧几里得，我们来看它与乘法逆元的关系。

  • 逆元：\(a\) 关于 模\(b\) 的逆元 整数\(d\) 满足 \(ad \bmod b \equiv 1\)
  • 扩展欧几里得：求方程 \(ax + by = GCD(a, b)\) 的一组可行解

逆元的\(ad \bmod b \equiv 1\)，等价于 \(ad-kb=1\)，其中\(k\)为未知整数。

设 \(d\) 为 \(x\)，\(-k\) 为 \(y\)，则\(ad-kb=1\)转换为 \(ax+by=1\)，

求出 \(x\) 就得到了 \(a\) 关于模 \(b\) 的逆元。

1
2
3
4
5
6
int ExGcdInv(int a, int b)
{
    int x, y;
    ExGCD(a, b, x, y);
    return x;
}

时间复杂度：大约O(logn)（斐波那契复杂度）。

适用范围：存在逆元即可求，适用于个数不多但模数b很大的时候，最常用、安全的求逆元方式。

方法2：费马小定理

除了扩展欧几里得，还有另一个方法可以求逆元。

费马小定理：对于整数 \(a\) 与质数 \(b\) ，若 \(a\) 与 \(b\) 互质，则有：

\[ a^{b − 1} \bmod b \equiv 1 \]

快速幂取模

快速幂取模大家应该也学过了，这里复习一下：

求x ^ n % MOD， n 很大时需要用折半的思想。如下所示求2^15：

1
2
3
4
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
4   4   4   4   4   4   4
16      16      16
256

可以看到，两两结合的时候，如果数字个数是奇数就会有“零头”，把零头存入ret，最终结果就是256*2*4*16。

1
2
3
4
5
6
7
8
9
10
11
int PowMod(int a, int n, int mod)
{
    int ret = 1;
    while(n)
    {
        if(n & 1) ret = ret * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return ret;
}

费马小定理求逆元

上文费马小定理的式子等价于

\[ a \cdot a^{b-2} \bmod b \equiv 1 \]

显然 \(a^{b-2}\) 就是 \(a\) 模 \(b\) 的逆元。

求逆元，就用 b-2 和 b 代替 快速幂取模中的 n 和 mod：

1
2
3
4
int FermatInv(int a, int b)
{
    return PowMod(a, b - 2, b);
}

时间复杂度：大约O(log b)。

适用范围：一般在模数 b 是质数的时候。

方法3：递归/递推求逆元

求 \(a\) 在模 \(b\) 时的逆元（如果 \(a > b\)，先将 \(a\) 取模 \(a:=a \bmod b\) 再求逆元，其中模数 \(b\) 是质数。

设 $k = b / a \(，\)r = b i$，则有

\(b = ak + r\)

=>

\(ak + r \equiv 0, (\bmod b)\)

=>

\(kr^{-1} + a^{-1} \equiv 0, (\bmod b)\)

=>

\(a^{-1} \equiv -kr^{-1}, (\bmod b)\)

\(a^{-1}\) 和 \(r^{-1}\) 可分别由其逆元 \(inv(a)\)和\(inv(r)\)代替，等式依然成立：

\(inv(a) \equiv -k \cdot inv(r), (\bmod b)\)

将 \(k\) 和 \(r\) 由 \(a\)、\(b\) 表达代入得：

\(inv(a) \equiv - \lfloor b / a \rfloor \cdot inv(b \bmod a), (\bmod b)\)

从而我们得到了由 \(inv(b \bmod a)\) 推出 \(inv(a)\) 的递推关系。

用递归方式计算的话，由于 \(b\) 是质数，\(b \bmod (b \bmod (b \bmod ...a))\) 总会到 \(1\)， \(inv(1) \equiv 1\)，从而达到递归终点。

1
2
3
4
5
LL Inv(LL a, LL b)
{
    if(a == 1) return 1;
    return (b - b / a) * Inv(b % a, b) % b;
}

时间复杂度：大约O(log b)

适用范围： \(b\) 一定要为质数，此方法代码简单，但使用需谨慎。

逆元打表

基于递推方法，就可以打表

1
2
3
4
5
6
7
int invList[mod + 10];
void GetInv(int mod)
{
    invList[1] = 1;
    for(int i = 2; i < mod; i ++)
        invList[i] = 1LL * (mod - mod / i) * invList[mod % i] % mod;
}

时间复杂度：显而易见。

适用场景：频繁调用不同数的逆元。