求逆元的四大基本方法
前几天在u裙见到有人说求乘法逆元四大方法,网上讲的全面的又很少,就来写一下(
乘法逆元:关于x的方程ax≡1(mod p)的解称为a模p的乘法逆元,记做a\(^{-1}\)
1、快速幂
仅限p是质数
a^(p-1)≡1(mod p),则a*a^(p-2)≡1(mod p),即a的逆元为a^(p-2)
代码
const int P=998244353;
int qpow(ll a,ll p){
ll r=1;
for(;p;p>>=1,a=a*a%P) if(p&1) r=r*a%P;
return r;
}
int inv(ll a){
return qpow(a,P-2);
}
2、欧几里得
要求a与p互质
考虑求a的逆元就是解一个方程ax≡1(mod b),即ax+by=1
这就很好做了,注意exgcd之后要把x搞成正的
代码
void exgcd(int a,int b,int &x,int &y){
if(!b){
x=1,y=0;
return;
}
exgcd(b,a%b,y,x);
y-=(a/b)*x;
}
int inv(int a,int p){
int x,y;
exgcd(a,p,x,y);
x=(x%b+b)%b;
}
3、递推
这个做法应该人人都会?
\(设ax+b=P,其中1<x<P,0\leq b<x\)
\(则ax+b \equiv 0(\bmod P)\)
同时乘上\(x^{-1}b^{-1}\)得到
\(ab^{-1}+x^{-1} \equiv 0(\bmod P)\)
\(x^{-1} \equiv -ab^{-1}(\bmod P)\)
代码
inv[1]=1;
for(int i=2;i<P;i++) inv[i]=((P-P/i)*inv[P%i])%P;
顺别提一嘴,可以用这个公式递归的求一个数的逆元,复杂度好像上界好像是O(n^1/3),期望O(logn)
参见OI-wiki 知乎相关讨论
手动分鸽-----------------------------------------------------------------------------------------------------------------------------------------
还有一种递推阶乘的逆元,大概就是预处理一下阶乘f,随便找个方法求出inv(f[n]),
然后for(int i=n-1;i>=1;i--) inv[i]=inv[i+1]*(i+1)%P;
inv[i]表示f[i]的逆元
手动分鸽-----------------------------------------------------------------------------------------------------------------------------------------
小拓展:O(n)求a1、a2、...、an的逆元
记\(f_i=a_1*a_2*...*a_i,g_i为f_i的逆元\)
用类似阶乘逆元的方式可以求出g_i
最后\(a_i的逆元就是g_i*f_{i-1}\)
4、牛顿迭代法
把求逆元看做求方程(ax-1)/x=0的解
根据牛顿迭代法列出柿子
就是说(注意,这里是对于实数)
\(
f(x)=\frac{ax-1}{x}\\
则f'(x)=\frac{1}{x^2}\\
x_{n+1}=x_n-\frac{f(x)}{f'(x)}=x_n-\frac{(ax_n-1)}{x_n}*x_n^2=x_n-x_n(ax_n-1)=x_n(2-ax_n)
\)
然后根据这篇文章[1]里的东西
懒得看的话大概意思就是说你有了
\(
ax_n \bmod 2^k = 1\)之后
可以得到
\(ax_{n+1} \bmod 2^{2k} = 1
\)
代码
ull inv(ull x){
ull r=1;
r*=2-r*x,r*=2-r*x,r*=2-r*x,r*=2-r*x,r*=2-r*x,r*=2-r*x;
return r;
}
这里的初值可以改变,例如令r=x是x mod 8的逆元,令r=x^2+x-1是x mod 16的逆元,而r=1就是x mod 2的逆元。最猛的可能是r=3x xor 2是x mod 32的逆元,具体参见上面的链接
最后的最后:
牛顿迭代法的推广可以直接去看这篇论文[2]
参考: