乘法逆元

数学基础

在模运算中，加减乘法都可以简单地进行转换，但是除法却不能这么做。例如 $(12÷2)mod7$，直接计算其值为 $6$，但是如果我们自作聪明，使用类似加减乘法的方式进行了转换 $(12mod7)÷2$，我们会发现结果变成了 $2.5$，这显然是不正确的。

在普通的计算中，除以一个数等于乘上这个数的倒数（即逆元），即 $a÷b=a\cdot \frac{1}{b}=a\cdot b^{-1}$。在模运算中这样做依然是可行的，$b^{-1}$ 就是 $b$ 的逆元，但是 $b^{-1}$ 的值不再是简单的求倒数，而是需要满足 $b\cdot b^{-1}modp=1$，记为 $b\cdot b^{-1}\equiv 1(modp)$。在之前的例子中，$2$ 在 $mod7$ 下的逆元是 $4$，所以 $(12÷2)mod7=(12\times 4)mod7=6$。

不同的 $p$ 下，$b^{-1}$ 的值不一定相同！

至此，如何在除法中引入模运算的问题被转化为了求取逆元的问题，而这一问题有非常多种解法，以下围绕公式 $a\cdot b\equiv 1(modp)$ 展开，$a$ 和 $b$ 互为 $modp$ 意义下的逆元。

费马小定理

费马小定理指出，当 $p$ 为素数且 $a$ 和 $p$ 互素时，有以下等式成立 $a^{p-1}\equiv 1(modp)$。

两边同时乘上 $a^{-1}$ 得 $a^{p-2}\equiv a^{-1}(modp)$，因此 $a^{-1}=a^{p-2}modp$，只需要使用快速幂取模算法即可迅速求出 $a$ 在 $modp$ 下的逆元。

费马小定理要求 $p$ 为素数，且 $a$ 和 $p$ 互素！

欧拉定理

欧拉函数 $\varphi (n)$ 表示 $[1,n]$ 中与 $n$ 互素的整数的个数，其定义为 $\varphi (n)=\sum _{i=1}^n[gcd(i,n)=1]$。显然当 $n$ 为素数时，$\varphi (n)=n-1$。

记 $k$ 为 $n$ 的质因子的个数，$facto{r}_n[i]$ 为 $n$ 的第 $i$ 个质因子，$i=1,2,\dots ,k$。这里直接给出关于欧拉函数的两个重要结论：

1. $\varphi (n)=n\cdot \prod _{i=1}^k(1-\frac{1}{facto{r}_n[i]})$
2. 当 $gcd(m_1,m_2)=1$ 时，$\varphi (m_1\cdot m_2)=\varphi (m_1)\cdot \varphi (m_2)$

给出欧拉定理，当 $gcd(a,p)=1$ 时， $a^{\varphi (p)}\equiv 1(modp)$。特别地，当 $p$ 为素数时，$a^{\varphi (p)}\equiv a^{p-1}\equiv 1(modp)$，即为费马小定理。

欧拉定理证明

1. 取 $modp$ 的缩系 $a_1,a_2,\dots ,a_{\varphi (p)}$，则 $a\cdot a_1,a\cdot a_2,\dots ,a\cdot a_{\varphi (p)}$ 也是 $modp$ 的缩系
2. 构造等式 $\prod _{i=1}^{\varphi (p)}{a}_i\equiv \prod _{i=1}^{\varphi (p)}a\cdot a_i\equiv a^{\varphi (p)}\cdot \prod _{i=1}^{\varphi (p)}{a}_i(modp)$
3. 等式两边同时除去 $\prod _{i=1}^{\varphi (p)}{a}_i$ ，即可得 $1\equiv a^{\varphi (p)}(modp)$

欧拉定理要求 $a$ 和 $p$ 互素。

费马小定理是欧拉定理在 $p$ 为素数时的特殊情况。

扩展欧几里得

根据裴蜀定理：若 $gcd(a,b)=d$ 且 $a,b\in \mathbb{Z}$，则对于任意 $x,y\in \mathbb{Z}$，都有 $(a\cdot x+b\cdot y)|d$。特别地，一定存在 $x,y\in \mathbb{Z}$，使得 $a\cdot x+b\cdot y=d$ 成立。

推导拓展欧几里

1. 假设 $a>b$
2. 显然当 $b=0$ 时
   1. $gcd(a,b)=a$
   2. 此时 $x=1,y=0$
3. 当 $a\cdot b\ne 0$ 时
   1. 设 $\{\begin{array}{ll}gcd(a,b)& =a\cdot x_1+b\cdot y_1\\ gcd(b,amodb)& =b\cdot x_2+(amodb)\cdot y_2\end{array}$
   2. 根据朴素的欧几里得原理有 $gcd(a,b)=gcd(b,amodb)$
   3. 则 $a\cdot x_1+b\cdot y_1=b\cdot x_2+(amodb)\cdot y_2=b\cdot x_2+(a-\lfloor a/b\rfloor \cdot b)\cdot y_2$
   4. 整理得 $a\cdot x_1+b\cdot y_1=a\cdot y_2+b\cdot (x_2-\lfloor a/b\rfloor \cdot y_2)$
   5. 可得 $x_1=y_2,y_1=x_2-\lfloor a/b\rfloor \cdot y_2$
4. 综上，可得 $x,y$ 的递推式并求解

公式 $a\cdot b\equiv 1(modp)$ 等价于 $a\cdot b+p\cdot k=1$，根据上述推理，若满足 $gcd(a,p)=1$，则可求解 $b,k$，其中 $b$ 即为 $a$ 在 $modp$ 下的逆元。

拓展欧几里得要求 $a$ 和 $p$ 互素！

线性递推

若 $p$ 为素数，线性递推法可以在 $O(p)$ 的时间内求出 $[1,p-1]$ 中每个整数在 $modp$ 下的逆元。

记 $inv[i]$ 为 $i$ 在 $modp$ 下的逆元，则有如下递推式：

$$inv[i]=\{\begin{array}{ll}1,& i=1\\ -\lfloor p/i\rfloor \cdot inv[pmodi],& i=2,3,\dots ,p-1\end{array}$$

递推式证明

1. 令 $a\cdot i+b=p$，其中 $a=\lfloor p/i\rfloor ,b=pmodi$
2. 有 $b\cdot inv[b]\equiv (p-a\cdot i)\cdot inv[b]\equiv (p\cdot inv[b]-a\cdot i\cdot inv[b])\equiv 1(modp)$
3. 因为 $p\cdot inv[i]\equiv 0(modp)$，有 $-a\cdot i\cdot inv[b]\equiv 1(modp)$
4. 代入 $a,b$ 得 $-\lfloor p/i\rfloor \cdot inv[pmodi]\cdot i\equiv 1\equiv i\cdot inv[i](modp)$
5. 所以 $inv[i]=-\lfloor p/i\rfloor \cdot inv[pmodi]$
6. 因为 $i>pmodi\ge 1$ 且 $inv[1]=1$，所以在知晓前 $k$ 项时，可以在常数时间内计算出第 $k+1$ 项，递推式成立

线性递推要求 $p$ 为素数，且 $a$ 和 $p$ 互素！

代码

幂运算的结果可能很大，应当注意数据类型以避免溢出！

费马小定理

Typescript

export function modInverseFermat(a: bigint, p: bigint): bigint {
  return quickMod(a, p - 2n, p);
}

Python

def mod_inverse_fermat(a: int, p: int) -> int:
    return quick_mod(a, p - 2, p)

ts#quickMod

/**
 * number类型支持的最大安全整数仅有2^53-1，大约为9e15
 * 在幂运算中非常容易溢出，建议使用更安全的bigint
 */
export function quickMod(x: bigint, p: bigint, mod: bigint): bigint {
  let res: bigint = 1n;
  while (p) {
    if (p & 1n) res = (res * x) % mod;
    x = (x * x) % mod;
    p >>= 1n;
  }
  return res;
}

py#quick_mod

def quick_mod(x: int, p: int, mod: int) -> int:
    res: int = 1
    while p:
        if p & 1:
            res = (res * x) % mod
        x = (x * x) % mod
        p >>= 1
    return res

欧拉定理

Typescript

export function eulerPhi(n: bigint) {
  let res: bigint = n;
  for (let i = 2n; i * i <= n; i++) {
    if (n % i === 0n) {
      res = (res / i) * (i - 1n);
      while (n % i === 0n) n /= i;
    }
  }
  if (n > 1n) res = (res / n) * (n - 1n);
  return res;
}

Python

def euler_phi(n: int) -> int:
    res = n
    for i in range(2, math.floor(math.sqrt(n)) + 1):
        if n % i == 0:
            res = (res // i) * (i - 1)
            while n % i == 0:
                n //= i
    if res > 1:
        res = (res // n) * (n - 1)
    return res

拓展欧几里得

Typescript

export function modInverseExgcd(a: number, p: number): number {
  const [gcd, x, y] = exgcd(a, p);
  if (gcd !== 1) return -1; // a和p不互质，无法求出逆元
  return (x + p) % p; // x可能为负数，需要加上p确保为正数
}

Python

def mod_inverse_exgcd(a: int, p: int) -> int:
    gcd, x, y = exgcd(a, p)
    if gcd != 1:
        return -1 # a和p不互质，无法求出逆元
    return (x + p) % p # x可能为负数，需要加上p确保为正数

ts#exgcd

export function exgcd(a: number, b: number): [number, number, number] {
  if (!b) return [a, 1, 0];
  const [gcd, x, y] = exgcd(b, a % b);
  return [gcd, y, x - Math.floor(a / b) * y];
}

py#exgcd

def exgcd(a: int, b: int) -> int:
    if b == 0:
        return a, 1, 0
    gcd, x, y = exgcd(b, a % b)
    return gcd, y, x - a // b * y

线性递推

Typescript

export function inverseLinear(p: number): number[] {
  const inv: number[] = new Array<number>(p).fill(0);
  inv[1] = 1;
  for (let i = 2; i <= p - 1; i++)
    inv[i] = ((p - Math.floor(p / i)) * inv[p % i]) % p;
  return inv;
}

const invCache: Map<number, number[]> = new Map();

export function inverseLazy(a: number, p: number): number {
  if (!invCache.has(p)) invCache.set(p, [0, 1]); // 初始化
  const array = invCache.get(p);
  if (array[a] === void 0)
    array[a] = ((p - Math.floor(p / a)) * inverseLazy(p % a, p)) % p;
  return array[a];
}

Python

def inverse_liner(p: int) -> List[int]:
    inv: List[int] = []
    inv[1] = 1
    for i in range(2, p):
        inv[i] = (p - p // i) * inv[p % i] % p
    return inv


@lru_cache(maxsize=None)
def inverse_lazy(a: int, p: int) -> int:
    if a == 1:
        return 1
    return (p - p // a) * inverse_lazy(p % a, p) % p