乘法逆元及其 python 求解法

定义

  乘法逆元是群论中的一个概念，它也在各类算法题中重现，是最好需要掌握的一项技巧。在同余中，乘法逆元的简单定义是，满足下式的 $x$ 称为 $a$ 模 $n$ 的乘法逆元。
$$ax\equiv 1(\mathrm{mod}n)\text{\hspace{1em}(1)}$$ 然而这个乘法逆元 $x$ 还不一定是存在的。$a$ 模 $n$ 的乘法逆元存在的充要条件是，$a$ 和 $n$ 互质。
  这个乘法逆元的作用就是，求 $y/a(\mathrm{mod}n)$ 可以转化为求 $yx(\mathrm{mod}n)$。

举一个实用的例子，求组合数 $C(n,k)=\frac{n!}{k!(n-k)!}$ 对一个很大的质数 $p$ 的模。显然，如果只是一个 $(n,k)$ 数据，暴力求一下也就行了。但如果给你很多很多 $(n,k)$ 的数据，时间复杂度就不可接受了。此时乘法逆元的作用就表现出来：
假设我们预先知道 $n$ 的上界，那我们就可以用 $O(n)$ 的复杂度构建fact数组，其中fact[i]的值是 $i!\mathrm{mod}p$。另外，再用 $O(n\log p)$ 的复杂度构建invs数组，其中invs[i]的值是 $i!$ 模 $p$ 的乘法逆元。
在经过 $O(n\log p)$ 复杂度的预处理之后，求解组合数就变成了 $O(1)$ 复杂度的事情。即C(n,k)=fact[n] * invs[k] * invs[n - k] % p。

gcd 算法

  先介绍一下求解逆元的前置知识。

辗转相除法求最大公因数

  相信大家在学习一门编程语言的时候，接触过求解最大公因数的算法。最出名的莫过于辗转相除法。记 $a,b$ 的最大公因数为 $\gcd (a,b)$，那么 $\gcd (a,b)=\gcd (b,a\mathrm{mod}b)$。
  当时你大概没有考虑过这个算法的时间复杂度是多少。实际上，使用辗转相除法求解 $a,b$ 的最大公因数，时间复杂度是 $O(\log (\max (a,b)))$。

拓展 gcd 算法

  先要用到拓展 gcd算法。这个算法的目标是，对于任意的正整数 $a,b$ 找到 $x,y$ 满足：
$$ax+by=\gcd (a,b)\text{\hspace{1em}(2)}$$   注意到当 $b=0$ 时，上式化为 $ax+by=b$，那么就找到了 $x=0,y=1$。
  而当 $b\ne 0$ 时，假设有一个 $\hat{x},\hat{y}$ 满足 $$b\hat{x}+(a\mathrm{mod}b)\hat{y}=\gcd (b,a\mathrm{mod}b)\text{\hspace{1em}(3)}$$ 根据辗转相除法我们知道 $(2)(3)$ 两式的右边是相等的，因此得到 $$ax+by=b\hat{x}+(a\mathrm{mod}b)\hat{y}\text{\hspace{1em}(3)}$$ 其中，令 $a=kb+t(k,t\in \mathbb{Z},0\le t<b)$，那么 $(3)$ 又可以化为 $$kbx+tx+by=b\hat{x}+t\hat{y}\text{\hspace{1em}(4)}$$ 稍微整理一下，可得 $$b(kx+y-\hat{x})+t(x-\hat{y})=0\text{\hspace{1em}(5)}$$   因此，假设我们已经求得了 $\hat{x},\hat{y}$，那么 $x=\hat{y}$ 和 $y=\hat{x}-kx=\hat{x}-\hat{y}(a//b)$ 就是我们想要一组 $(x,y)$ 的值。其中 $a//b=\lfloor a/b\rfloor (a,b>0)$。
  由拓展 gcd 算法中，求解 $(x,y)$ 和求 $\gcd (a,b)$ 是同步的，辗转相除法求解 $\gcd (a,b)$ 的终止条件也是 $b=0$。因此拓展 gcd 算法的时间复杂度也应该是 $O(\log (\max (a,b)))$，这个 gcd 在拓展算法中起到一个引路的作用。

求解模素数的乘法逆元

  有了上面的基础，就可以求解模素数的乘法逆元了。假设 $n$ 是一个正素数，$a>0$ 且 $a\ne n$，那我们可以在 $O(\log (\min (a,n)))$ 的时间复杂度内求得满足 $ax+ny=\gcd (a,n)=1$ 的 $(x,y)$。式子两边同时对 $n$ 去模，得到 $ax\equiv 1(\mathrm{mod}n)$。因此，通过拓展 gcd 算法求得的 $x$ 就是我们想要的乘法逆元。
  附带 python 代码一份，用于求 $b$ 模 $k$ 的乘法逆元。

def extended_gcd(a, b):
    if a == 0:
        return b, 0, 1
    gcd, x1, y1 = extended_gcd(b % a, a)
    x = y1 - (b // a) * x1
    y = x1
    return gcd, x, y

def mod_inverse(b, k):
    gcd, x, _ = extended_gcd(b, k)
    if gcd != 1:
        raise ValueError(f"No modular inverse for {b} mod {k}")
    return x % k