// Extended Euclidean algorithm.void ex_gcd(int a, int b, int& x, int& y) { if (!b) { x = 1; y = 0; } else { ex_gcd(b, a % b, y, x); y -= a / b * x; }}// Returns the modular inverse of a modulo m.// Assumes that gcd(a, m) = 1, so the inverse exists.int inverse(int a, int m) { int x, y; ex_gcd(a, m, x, y); return (x % m + m) % m;}
Python
# Extended Euclidean algorithm.def ex_gcd(a, b): if b == 0: return 1, 0 else: x1, y1 = ex_gcd(b, a % b) x = y1 y = x1 - (a // b) * y1 return x, y# Returns the modular inverse of a modulo m.# Assumes that gcd(a, m) = 1, so the inverse exists.def inverse(a, m): x, y = ex_gcd(a, m) return (x % m + m) % m
// Binary exponentiation.int pow(int a, int b, int m) { long long res = 1, po = a; for (; b; b >>= 1) { if (b & 1) res = res * po % m; po = po * po % m; } return res;}// Returns the modular inverse of a prime modulo p.int inverse(int a, int p) { return pow(a, p - 2, p); }
Python
# Returns the modular inverse of a prime modulo p.# Use built-in pow function.def inverse(a, p): return pow(a, p - 2, p)
当然,理论上,这一方法可以利用 欧拉定理 推广到一般的模数 m 的情形,即利用 aφ(m)−1modm 计算逆元。但是,单次求解 欧拉函数φ(m) 并不容易,因此该算法在一般情况下效率不高。
多个逆元的求法
有些场景下,需要快速处理出多个整数 a1,a2,⋯,an 在模 m 意义下的逆元。此时,逐个求解逆元,总共需要 O(nlogm) 的时间。实际上,如果将它们统一处理,就可以在 O(n+logm) 的时间内求出所有整数的逆元。
考虑序列 {ai} 的前缀积:
S0=1,Si=aiSi−1,i=1,2,⋯,n.
只要每个 ai 都与 m 互素,它们的乘积 Sn 就与 m 互素。因此,可以通过前文所述算法求出 Sn−1modm 的值。因为乘积的逆元就是逆元的乘积,所以,从 Sn−1 出发,反向遍历序列就能求出每个 Si 的逆元:
Si−1−1=aiSi−1modm,i=n,n−1,⋯,1.
由此,单个 ai 的逆元可以通过下式计算:
ai−1=Si−1Si−1modm,i=1,2,⋯,n.
参考实现如下:
参考实现
[list2tab]
C++
// Returns the modular inverses for each x in a modulo m.// Assume x mod m exists for each x in a.std::vector<int> batch_inverse(const std::vector<int>& a, int m) { int n = a.size(); std::vector<int> prod(n); long long s = 1; for (int i = 0; i < n; ++i) { // prod[i] = product of a[0...i-1]; prod[0] = 1. prod[i] = s; s = s * a[i] % m; } // s = product of all elements in a. s = inverse(s, m); std::vector<int> res(n); for (int i = n - 1; i >= 0; --i) { res[i] = s * prod[i] % m; s = s * a[i] % m; } return res;}
Python
# Returns the modular inverses for each x in a modulo m.# Assume x mod m exists for each x in a.def batch_inverse(a, m): n = len(a) prod = [0] * n s = 1 for i in range(n): # prod[i] = product of a[0...i-1]; prod[0] = 1. prod[i] = s s = s * a[i] % m # s = product of all elements in a. s = inverse(s, m) res = [0] * n for i in reversed(range(n)): res[i] = s * prod[i] % m s = s * a[i] % m return res
算法中,只求了一次单个元素的逆元,因此总的时间复杂度是 O(n+logm) 的。
线性时间预处理逆元
如果要预处理前 n 个正整数在素数模 p 下的逆元,还可以通过本节将要讨论的递推关系在 O(n) 时间内计算。这一方法常用于组合数计算中前 n 个正整数的阶乘的倒数的预处理。
对于 1<i<p 的正整数 i,考察带余除法:
p=⌊ip⌋i+(pmodi).
将该等式对素数 p 取模,就得到
0≡⌊ip⌋i+(pmodi)(modp).
将等式两边同时乘以 i−1(pmodi)−1 就得到
i−1≡−⌊ip⌋(pmodi)−1(modp).
这就是用于线性时间递推求逆元的公式。由于 pmodi<i,这一公式将求解 i−1modp 的问题转化为规模更小的问题 (pmodi)−1modp。因此,从 1−1modp=1 开始,对每个 i 顺次应用该公式,就可以在 O(n) 时间内获得前 n 个整数的逆元。
参考实现如下:
参考实现
[list2tab]
C++
// Precomputes modular inverses of all integers from 1 to n modulo prime p.std::vector<int> precompute_inverses(int n, int p) { std::vector<int> res(n + 1); res[1] = 1; for (int i = 2; i <= n; ++i) { res[i] = (long long)(p - p / i) * res[p % i] % p; } return res;}
Python
# Precomputes modular inverses of all integers from 1 to n modulo prime p.def precompute_inverses(n, p): res = [0] * (n + 1) res[1] = 1 for i in range(2, n + 1): res[i] = (p - p // i) * res[p % i] % p return res
这一算法只适用于模数是素数的情形。对于模数 m 不是素数的情形,无法保证递推公式中得到的 mmodi 仍然与 m 互素,因而递推所需要的 (mmodi)−1 可能并不存在。一个这样的例子是 m=8,i=3。此时,mmodi=2,不存在模 m 的逆元。
另外,得到该递推公式后,一种自然的想法是直接递归求解任意一个数 a 的逆元。每次递归时,都利用递推公式将它转化为更小的余数 pmoda 的逆元,直到余数变为 1 时停止。目前尚不清楚这样做的复杂度 1,因此,推荐使用前文所述的常规方法求解。