乘法逆元
含义
首先,数学上的乘法逆元就是指直观的倒数,即 a a a 的逆元是 1 a \displaystyle\frac{1}{a} a1 ,也即与 a a a 相乘得 1 1 1 的数。 a x = 1 ax=1 ax=1,则 x x x 是 a a a 的乘法逆元。
如果到度娘上搜索逆元,会看到如下定义:(这里逆元素是逆元的全称)
逆元素是指一个可以取消另一给定元素运算的元素 。
这个定义似乎不太那么直观……
我们来举个例子吧,先再实数范围举例,由小学知识可知,如果一个代数式 F F F 乘一个数 a a a 后,再乘它的倒数 1 a \frac 1 a a1 ,相当于没有乘 a a a (这里不考虑 0 0 0 的情况),换句话说,我们乘 1 a \frac 1 a a1 后,取消了代数式 F F F 乘 a a a 后值增大的影响。
不难发现这符合逆元的定义,故我们可以说一个数和其倒数互为乘法逆元。
而算法这里这里我们讨论关于取模运算的乘法逆元,即对于整数
a
a
a ,与
a
a
a 互质的数
b
b
b 作为模数,当整数
x
x
x 满足
a
x
m
o
d
b
≡
1
ax\mod\,b≡1
axmodb≡1 时,称
x
x
x 为
a
a
a 关于模
b
b
b 的逆元,代码表示是 a * x % b == 1
。
我们通常将 a a a 的逆元记作 a − 1 a^{-1} a−1 。
在算法竞赛中,经常会遇到求解数据很大,则输出模 1 0 9 + 7 10^9+7 109+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
。
故而我们需要一个算法求 除数 的 取模逆元 ,从而在四则运算取模的任务中,用逆元将除法转为乘法。
求法
1.费马小定理
前置芝士:欧拉定理:若 a , n a,n a,n互质,则 a ϕ ( n ) ≡ 1 ( m o d n ) a^{\phi(n)}\equiv 1 (\mod n) aϕ(n)≡1(modn)。
费马小定理:当有两数 a a a, p p p 满足 g c d ( a , p ) = 1 gcd ( a , p ) = 1 gcd(a,p)=1, 并且 p p p 是质数时,则有
a p − 1 ≡ 1 ( m o d p ) \begin{aligned} a^{p-1}\equiv 1\,\;(\mod p\;) \end{aligned} ap−1≡1(modp)
将该等式与
a
x
≡
1
(
m
o
d
p
)
ax\equiv 1(\mod p)
ax≡1(modp) 对比可以惊奇地发现:
x
=
a
p
−
2
\begin{aligned} x=a^{p-2} \end{aligned}
x=ap−2
逆元求出来之后,就可以用 快速幂 求出 a p − 2 a^{p-2} ap−2 的值了;
时间复杂度:大约 O ( l o g p ) O(log p) O(logp) 。
这里没啥好贴的代码了,贴一个快速幂吧:
//计算a^k%p
int qmi(int a,int k,int p){
int res=1;
while(k){
if(k&1) res=(LL)res*a%p;
k>>=1;
a=(LL)a*a%p;
}
return res;
}
2.扩展欧几里得
前置芝士:裴蜀定理:对于任意的正整数 a , b a,b a,b 一定存在正整数 x , y x,y x,y 使得 a x + b y = g c d ( a , b ) ; ax+by=gcd(a,b); ax+by=gcd(a,b); (并且也是能凑出来的最小正整数)
由定义可知: a x ≡ 1 ( m o d p ) ax\equiv 1\,(\mod p) ax≡1(modp),这个式子等价于已知 a , p a,p a,p 求一个二元一次不定方程 a x = k p + 1 ax=kp+1 ax=kp+1,移一下项得: a x − k p = 1 ax−kp=1 ax−kp=1。这东西与 a x + b y = g c d ( a , b ) ax+by=gcd(a,b) ax+by=gcd(a,b) 不是神似?就是扩展欧几里得算法。
附一小段扩展欧几里得代码:
//返回gcd(a,b),计算x,y
int exgcd(int a,int b,int &x,int &y){
if(!b){
x=1,y=0;
return a;
}
int d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
相比于费马小定理,扩展欧几里得的使用限制更小,不再需要 p p p 是质数,只需要 a , p a,p a,p 互质即可。
时间复杂度:大约 O ( l o g n ) O(logn) O(logn) 。
3.递推计算阶乘的逆元
当我们要计算一大串连续的阶乘的逆元时,采用费马小定理或扩展欧几里得算法就有可能超时,所以我们必须采用一个更快的算法。
令
f
i
=
i
!
f_i=i!
fi=i! ,则可得:
i
n
v
(
f
i
)
≡
i
n
v
(
f
i
−
1
⋅
(
i
)
)
(
m
o
d
p
)
\begin{aligned} inv(f_{i})\equiv inv(f_{i-1}\cdot (i))(\mod p) \end{aligned}
inv(fi)≡inv(fi−1⋅(i))(modp)
我们将右式乘开,则有:
i
n
v
(
f
i
)
≡
i
n
v
(
f
i
−
1
)
⋅
i
−
1
(
m
o
d
p
)
\begin{aligned} inv(f_{i})\equiv inv(f_{i-1})\cdot i^{-1}(\mod p) \end{aligned}
inv(fi)≡inv(fi−1)⋅i−1(modp)
最后我们每一次只用处理i这一个数的逆元了。
下面是一小段预处理代码(提前打表):
//fact[i]表示i的阶乘,infact[i]表示i的阶乘的逆元。
fact[0]=infact[0]=1;
for(int i = 1;i < N;i ++ ){
fact[i]=(LL)fact[i - 1]*i%mod;
infact[i]=(LL)infact[i - 1]*qmi(i,mod - 2,mod)%mod;
}
4.递推计算连续的数的逆元 (线性求逆)
当我们要计算连续的数的逆元时,我们可以采用以下递推式:
i
n
v
(
i
)
=
(
p
−
⌊
p
i
⌋
)
×
i
n
v
(
p
m
o
d
i
)
m
o
d
p
\begin{aligned} inv(i)=(p-\lfloor\frac{p}{i}\rfloor)\times inv(p\mod i)\mod p \end{aligned}
inv(i)=(p−⌊ip⌋)×inv(pmodi)modp
证明:设
t
=
⌊
p
i
⌋
,
k
=
p
m
o
d
i
t=\lfloor\frac{p}{i}\rfloor,k=p\mod i
t=⌊ip⌋,k=pmodi ,那么显然有
t
×
i
+
k
≡
0
(
m
o
d
p
)
\begin{aligned}t\times i+k\equiv 0(\mod p)\end{aligned}
t×i+k≡0(modp)
变形可得
−
t
×
i
≡
k
(
m
o
d
p
)
\begin{aligned} -t\times i\equiv k(\mod p) \end{aligned}
−t×i≡k(modp)
两边同时除以
i
k
ik
ik 得到
−
t
×
1
k
≡
1
i
(
m
o
d
p
)
\begin{aligned} -t\times\frac{1}{k}\equiv\frac{1}{i}(\mod p) \end{aligned}
−t×k1≡i1(modp)
即
i
n
v
(
i
)
≡
−
t
×
i
n
v
(
k
)
(
m
o
d
p
)
\begin{aligned} inv(i)\equiv-t\times inv(k)(\mod p) \end{aligned}
inv(i)≡−t×inv(k)(modp)
所以有
i
n
v
(
i
)
=
(
p
−
⌊
p
i
⌋
)
×
i
n
v
(
p
m
o
d
i
)
m
o
d
p
\begin{aligned} inv(i)=(p-\lfloor\frac{p}{i}\rfloor)\times inv(p\mod i)\mod p \end{aligned}
inv(i)=(p−⌊ip⌋)×inv(pmodi)modp
参考代码:
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
const int maxn = 3e6 + 10;
long long infact[maxn], n, p;
int main()
{
scanf("%lld%lld", &n, &p);
infact[0] = 0, infact[1] = 1;//初始化1的逆元为1
printf("%lld\n", infact[1]);
for (int i = 2; i <= n; ++i)
{
infact[i] = (p - p / i) * infact[p % i] % p;//上文中的式子
printf("%lld\n", infact[i]);
}
return 0;
}