乘法逆元

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

乘法逆元

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

这里我们讨论关于取模运算的乘法逆元，即对于整数 $a$，与 $a$ 互质的数 $b$ 作为模数，当整数 $x$ 满足 $axmodb\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/amodb$$

对于 $a$ 的逆元 $d$，能够满足

$$m\cdot dmodb=m/amodb=c$$

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

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

方法1：扩展欧几里得

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

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

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

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

设两个式子

$$ax+by=GCD(a,b)$$

$$b{x}^{\prime }+(amodb)y^{\prime }=GCD(b,amodb)$$

由欧几里得算法知 $GCD(a,b)=GCD(b,amodb)$

所以

$$ax+by=b{x}^{\prime }+(amodb)y^{\prime }$$

而

$$amodb=a-kb$$

其中 $k=\lfloor a/b\rfloor$，$\lfloor \rfloor$表示向下取整。

所以有

$$ax+by=b{x}^{\prime }+(a-kb)y^{\prime }$$

展开移项得：

$$ax+by=a{y}^{\prime }+b(x^{\prime }-k{y}^{\prime })$$

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

$x=y^{\prime }$、$y=x^{\prime }-k{y}^{\prime }$ 来进一步求一个可行解。

递归地用 $x^{\prime }$、$y^{\prime }$ 表示“上一步”的$x$、$y$，就能递归地把问题转换成

$$b{x}^{\prime }+(amodb)y^{\prime }=GCD(b,amodb)$$

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

$$x=y^{\prime }$$

$$y=x^{\prime }-k{y}^{\prime }$$

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

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$ 满足 $admodb\equiv 1$
• 扩展欧几里得：求方程 $ax+by=GCD(a,b)$ 的一组可行解

逆元的$admodb\equiv 1$，等价于 $ad-kb=1$，其中$k$为未知整数。

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

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

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}modb\equiv 1$$

快速幂取模

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

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

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。

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}modb\equiv 1$$

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

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

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

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

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

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

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

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

$b=ak+r$

=>

$ak+r\equiv 0,(modb)$

=>

$k{r}^{-1}+a^{-1}\equiv 0,(modb)$

=>

$a^{-1}\equiv -k{r}^{-1},(modb)$

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

$inv(a)\equiv -k\cdot inv(r),(modb)$

将 $k$ 和 $r$ 由 $a$、$b$ 表达代入得：

$inv(a)\equiv -\lfloor b/a\rfloor \cdot inv(bmoda),(modb)$

从而我们得到了由 $inv(bmoda)$ 推出 $inv(a)$ 的递推关系。

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

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

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

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

逆元打表

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

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;
}

时间复杂度：显而易见。

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