学这个之前建议看看之前发的那篇关于快速幂和扩展gcd的博文:跳转
一.引出逆元的用途
在讲逆元是什么之前,先说一下我们遇到的错误。。
对于特别大的数的连乘,让我们对最终结果取模1e9 + 7,可能在连乘的时候就已经爆longlong了,我们应该怎么做呢?
根据同余定理可知,(a✖️b) % mod = ( (a % mod) ✖️ (b % mod)) % mod,因此,我们可以设longlong型数ans = 1,每次乘以一个较大数的时候,都取模一下1e9+7,避免爆longlong。代码如下:
const ll mod = 1e9 + 7;
ll ans = 1;
for (int i = 0; i < n; i++) {
ans = ans * a[i] % mod;
}
最终ans即为n个数连乘取模1e9 + 7的结果
若是不只是乘法,中间再夹杂着除法,我们再按照上面的方式ans = ans / a[i] % mod,当我们兴冲冲的上交后,一个冷冰冰的WA在我们脸上胡乱的拍~
这是为啥呢?
因为同余定理中不包含(a / b) % mod = ((a % mod) / (b % mod)) % mod,记住不包含!!!
那么我们该怎么办呢?好办,我们将除法改成乘不就完了嘛~
除以一个数不就是乘一个数的倒数嘛,一个倒数在取模mod时为多少呢?
这时候逆元就出现了!!
也就是说,a / b % mod相当于a ✖️ b关于mod的逆元 % mod
我们只需要求b关于mod的逆元就可以了~
二.逆元的求法
至于求逆元的方式有很多,比如利用费马小定理,扩展欧几里得,欧拉定理等等…我比较常用方法是费马小定理和扩展欧几里得
(1)利用费马小定理求逆元:
我们由于费马小定理可知:当p为素数时,有a ^ (p - 1) ≡ 1 (mod p)
两边同时除a,得a^(p - 2) ≡ inv(a) (mod p) 其中inv(a)是a关于p的逆元
因此: inv(a) = a ^ (p - 2) (mod p)
我们就可以通过使用快速幂求得a ^ (p - 2)取模p得到 a 关于p的逆元
代码如下:
//利用快速幂求a^(p - 2) mod p
ll pow_mod (ll a, ll b, ll p) {
ll ans = 1;
while (b) {
if (b & 1) ans = ans * a % p;
b >>= 1;
a = a * a % p;
}
return ans;
}
//费马小定理求a关于p的逆元
ll permat (ll a, ll p) {
return pow_mod (a, p - 2, p);
}
(2)利用扩展欧几里得求逆元
上一篇我们学了gcd的求法,它是根据欧几里得算法的核心 gcd (a, b) = gcd (b, a % b)递归得到的
那么,扩展欧几里得是什么呢?
我们可以利用已知的a,b,一定存在至少一组解x,y,使它们满足贝祖等式:a✖️x + b✖️y = gcd(a, b)
至于这个定理为什么成立,我也没学。。。Orz,我觉得会用就可以了,冷汗?
根据a✖️x + b✖️y = gcd(a, b),若a,b互质,那么此时a✖️x + b✖️y = gcd(a, b) = 1
我们对a✖️x + b✖️y = 1同时取余b,会得到a✖️x % b + b✖️y % b = 1 % b
=> a✖️x % b = 1 % b
=> a✖️x = 1(mod b)
大家看着是不是似成相识呢?,对!x就是a关于b的逆元
因此,我们就可以通过扩展欧几里得算法求解x即可。
关于扩展欧几里得算法:
我们设x✖️a + y✖️b = d (d为gcd(a, b)); (1)
那么x1✖️b + y1✖️(a % b) = d; (2)
也就是说x✖️a + y✖️b = x1✖️b + y1✖️(a % b); (3)
因为a % b = a - (a / b)✖️b; (4)
我们把(4)带入(3),得:
x✖️a + y✖️b = x1✖️b + y1✖️(a - (a / b)✖️b) (5)
化简得:x✖️a + y✖️b = y1✖️a + (x1 - (a / b)✖️y1)✖️b (6)
因此,我们可得x = y1, y = x1 - (a / b)✖️y1
代码实现如下:
//扩展欧几里得
ll extend_gcd (ll a, ll b, ll &x, ll & y) {
if (a == 0 && b== 0) return -1; //此时由于a,b无最大公约数,所以出错
if (b == 0) { //递归终止条件,若b == 0,那么gcd(a,b) == a,a * x + b * y = gcd (a, b) = 1,因此x = 1
x = 1;
y = 0;
return a; //返回a用于判断a是否为1,用以确定原a,b是否互质
}
ll d = extend_gcd (b, a % b, y, x); //由扩展欧几里得可得
y -= a / b * x;
return d;
}
//求逆元 ax = 1 (mod p)
ll mod_reverse (ll a, ll p) {
ll x, y;
ll d = extend_gcd(a, p, x, y); //把a带入a,p带入b,通过扩展欧几里得求x
if (d == 1) return (x % p + p) % p; //如果a == 1表明gcd(a,p)为1,也就是a,p互质,因此输出x
else return -1; //否则a,p不互质,返回-1
}
三.线性求逆元
求素数有线性筛,那么逆元呢?逆元当然不能落后啦,于是线性求逆元出现了Orz~
线性求逆元的公式是:inv(a) = (p - p / a)✖️inv(p % a) % p
证明过程:设x = p % a, y = p / a;
则有x + y✖️a = p
因此(x + y✖️a) % p = 0
x % p = (- y)✖️a % p
x✖️inv(a) % p = (-y) % p
inv(a) = (p - y)✖️inv(x) % p
inv(a) = (p - p / a)✖️inv(p % a) % p
这样,我们就可以利用之前的逆元求解后面的逆元了,代码如下:
void init() {
inv[0] = 1;
inv[1] = 1;
for (int i = 2; i <= N; i++) {
inv[i] = ((p - p / i) * inv[p % i]) % p;
}
}
四.线性求阶乘逆元
如果我们需要求0!到n!的逆元,对每个元素都求一遍会特别慢
前面说了,逆元就可一看做是求倒数
那么就有1 / (n+1)! × (n+1)=1/ n!
因此inv[n + 1]✖️(n + 1) = inv[n] (mod p)
代码如下:
ll fact[N + 5]; //存储阶乘
ll inv[N + 5]; //存储阶乘的逆元
fact[0] = 1;
for (int i = 1; i <= N; i++) {
fact[i] = fact[i - 1] * i % p; //线性求阶乘
}
inv[N] = mod_reverse(fact[N], p); //利用前面学的扩展欧几里得求fact[N]关于p的逆元
for (int i = N - 1; i >= 0; i--) {
inv[i] = inv[i + 1] * (i + 1) % mod; //求线性阶乘逆元
}
以上就是我对逆元的全部理解了,就酱紫~
如果有写的不对或者不全面的地方 可通过主页的联系方式进行指正,谢谢