乘法逆元
nyt_nyt2012
·
·
算法·理论
逆元是很多数论问题的基础，本文将介绍模意义下的乘法逆元，并讨论快速求出多个数逆元的方法。
概念
模意义下的乘法逆元主要用于计算 \frac{a}{b}\mod p 的值，a 的逆元记作 a^{-1}，且满足 a\times a^{-1}\equiv1\mod p。只有 \gcd(a,p)=1，即 a、p 互质时 a 才会有模 p 意义下的逆元。由于模意义下除法会出错，这时我们就需要引入逆元了，因为 \frac{a}{b}\equiv a\times b^{-1}\mod p，这样就把除法转换为了乘法，就可以计算了。
逆元的计算方式有很多，复杂度和适应场景也各不相同，让我们来看看吧！
逆元的求法
先简单说说单个数逆元的求法，不是本文重点，也比较简单，可以跳过。
:::info[单个数逆元的求法]
0. 暴力
你就说这是不是逆元的求法吧！
首先，a^{-1} 大于 p 是无用的，于是我们可以枚举 1\sim p-1，并验证其是不是真正的逆元。
int get_inv(int a,int p){
int inv=-1;
for(int i=1;i<p;i++)
if(a*i%p==1){
inv=i;
break;
}
return inv;
}
## 1. 快速幂
这个做法需要用到费马小定理，所以只有 $p$ 是质数时才能使用。
> 若 $p$ 为质数，$a$ 为正整数，且 $a$、$p$ 互质，则 $a^{p-1}\equiv1\mod p$。
由此我们可以得出：$a\times a^{p-2}\equiv1\mod p$，所以 $a^{p-2}\bmod p$ 就是 $a^{-1}$，这可以用快速幂算法在 $O(\log p)$ 的时间复杂度内完成一个数逆元的计算。
```cpp
int get_inv(int a,int p){
int inv=1,base=a,b=p-2;
while(b){
if(b&1)inv=inv*base%p;
base=base*base%p;
b>>=1;
}
return inv;
}
```
大部分比赛上，模数通常是质数，这种逆元求法常数小且好写好调，十分实用。
## 2. 扩展欧几里得
由于 $a\times a^{-1}\equiv1\mod p$，所以一定有整数 $y$ 使得 $a^{-1}\times a+y\times p=1=\gcd(a,p)$，这就是一个经典的线性同余方程，可以使用 exgcd 算法解决。这个方法的优点是可以处理任意模数，常数也很小，可以在 $O(\log\max(a,p))$ 的时间复杂度内完成一个数逆元的计算。
```cpp
void exgcd(int a,int b,int& x,int& y){
if(b==0){
x=1;y=0;
return;
}
exgcd(b,a%b,x,y);
int t=x;x=y;
y=t-a/b*y;
}
int get_inv(int a,int p){
int x=0,y=0;
exgcd(a,p,x,y);
return (x%p+p)%p;
}
```
:::
## 3. 线性求 $1\sim n$ 的逆元
我们默认 $n<p$，且 $1\sim n$ 的整数都与 $p$ 互质。
对于某个数 $i$，设 $p=i\times a+b$。我们就可以得到：
$$
i\times a+b\equiv0\mod p
$$
发现可以进行转化，两边同时乘以 $i^{-1}\times b^{-1}$ 会有奇效：
$$
a\times b^{-1}+i^{-1}\equiv0\mod p
$$
由于 $b<i,1^{-1}=1$，所以上式可以递推计算，时间复杂度为 $O(n)$。
```cpp
vector<int> get_inv(int n,int p){
vector<int>inv(n+1);
inv[1]=1;
for(int i=2;i<=n;i++)inv[i]=inv[p%i]*(p-p/i)%p;
return inv;
}
```
这是封装到一个函数里的写法，可能在[模板题](/problem/P3811)的代码更易读：
:::success[Code]
```cpp
#include<bits/stdc++.h>
#define int long long
using namespace std;
int inv[3003030],n,p;
signed main(){
ios::sync_with_stdio(0);
cin.tie(0);
cin>>n>>p;
inv[1]=1;
for(int i=2;i<=n;i++)inv[i]=inv[p%i]*(p-p/i)%p;
for(int i=1;i<=n;i++)cout<<inv[i]<<'\n';
return 0;
}
```
:::
## 4. 离线求任意 $n$ 个数的逆元
这个算法利用了逆元的可乘性，即 $a^{-1}\times b^{-1}\equiv (a\times b)^{-1}\mod p$。
首先，我们离线读入待求逆元的数 $a_1\sim a_n$，我们可以算出其前缀积和后缀积，并计算 $\Pi^n_{i=1}a_i$ 的逆元，由逆元的可乘性，就有：
$$
(\Pi^n_{i=1}a_i)^{-1}\equiv\Pi^n_{i=1}a_i^{-1}\mod p
$$
又因为 $a\times a^{-1}\equiv1\mod p$，所以任意一段连续的数的乘积的逆元都可以通过上式乘以某个前缀积和某个后缀积计算，计算方式如下：
$$
(\Pi^r_{i=l}a_i)^{-1}\equiv\Pi^n_{i=1}a_i^{-1}\times\Pi^{l-1}_{i=1}a_i\times\Pi^n_{i=r+1}a_i\mod p
$$
对于单个数的逆元，可以看做长度为 $1$ 的一个区间的乘积的逆元，可以同上计算，时间复杂度为 $O(n+\log p)$。
```cpp
vector<int> get_inv(int n,int a[],int p){
vector<int>inv(n+2),qz(n+2),hz(n+2);
qz[0]=hz[n+1]=1;
for(int i=1;i<=n;i++)qz[i]=qz[i-1]*a[i]%p;
for(int i=n;i>=1;i--)hz[i]=hz[i+1]*a[i]%p;
int invall=1,base=qz[n],b=p-2;
while(b){
if(b&1)invall=invall*base%p;
base=base*base%p;
b>>=1;
}
for(int i=1;i<=n;i++)inv[i]=invall*qz[i-1]%p*hz[i+1]%p;
return inv;
}
```
同样的，这是封装到一个函数里的写法，可能在[例题](/problem/P5431)的代码更易读：
:::success[Code]
```cpp
#include<bits/stdc++.h>
#include<sys/mman.h>
#include<sys/stat.h>
#include<fcntl.h>
#include<unistd.h>
using namespace std;
const int maxn=5e6+11;
int k,a[maxn],n,inv[maxn],qz[maxn],hz[maxn],p;
uint8_t* str;
inline int read(){
int k=0;uint8_t c=(*str++);
while(!isdigit(c))c=(*str++);
while(isdigit(c))k=k*10+(c&15),c=(*str++);
return k;
}
signed main(){
struct stat ss;
fstat(0,&ss);
str=(uint8_t*)mmap(NULL,ss.st_size,1,2,0,0);
n=read();p=read();k=read();
for(int i=1;i<=n;i++)a[i]=read();
qz[0]=hz[n+1]=1;
for(int i=1;i<=n;i++)qz[i]=(long long)qz[i-1]*a[i]%p;
for(int i=n;i>=1;i--)hz[i]=(long long)hz[i+1]*a[i]%p;
long long invall=1,base=qz[n],b=p-2;
while(b){
if(b&1)invall=invall*base%p;
base=base*base%p;
b>>=1;
}
for(int i=1;i<=n;i++)inv[i]=invall*qz[i-1]%p*hz[i+1]%p;
long long ans=0;base=1;
for(int i=1;i<=n;i++){
base=base*k%p;
ans=(ans+base*inv[i])%p;
}
cout<<ans;
return 0;
}
```
同时，上代码使用了快速幂的方式求所有数乘积的逆元，只能处理模数是质数的情况。如果需要在任意模数下使用该方法需要用扩欧求所有数乘积的逆元。
:::
## 5. 在线求逆元
重点来了！在线 $O(1)$ 逆元！
定义 Farey 序列 $F_n$ 为母不超过 $n$ 的真分数和 $\frac01$ 以及 $\frac11$ 组成的序列，例如，$F_5$ 如下：
$$
\bigg\{\frac01,\frac15,\frac14,\frac13,\frac25,\frac12,\frac35,\frac23,\frac34,\frac45,\frac11\bigg\}
$$
有一个引理：
> 对于任意整数 $m\ge2$ 和任意实数 $v\in[0,1]$，总能在 $F_{m-1}$ 阶 Farey 序列中找到一个分数 $\frac{x}{y}$，使得 $\lvert v-\frac{x}{y}\rvert\le\frac1{y\times m}$。
>
设 $v=\frac{a}{b}$，则有：
$$
\lvert\frac{a}{b}-\frac{x}{y}\rvert\le\frac1{y\times m}
$$
对不等式两边同时乘以 $b\times y$，就会得到：
$$
\lvert a\times y-b\times x\rvert\le\frac{b}{m}
$$
把 $b$ 换成模数 $p$，则 $a\times y\equiv k\mod p$，其中 $\lvert k\rvert\le\lfloor\frac{p}{m}\rfloor$，这里我们可以预处理出 $1\sim \lfloor \frac{p}{m}\rfloor$ 的逆元，之后通过其他方法规避负数问题。
问题来了，怎么找到这个 $\frac{x}{y}$ 呢？显然，$F_m$ 中相邻两项之差一定大于 $\frac1{m^2}$，故 $\lfloor\frac{m^2\times x}{y}\rfloor$ 一定互不相同，于是我们可以记录 $\lfloor\frac{m^2\times x}{y}\rfloor$ 中 $0\sim m^2$ 是否出现过，最后分别从前往后和从后往前记录前驱后继即可。
考虑上面这个鬼的时间复杂度，我们需要预处理出 $1\sim \lfloor \frac{p}{m}\rfloor$ 的逆元，这个可以做到 $O(\frac{p}{m})$；我们还需要预处理出 $\lfloor\frac{m^2\times x}{y}\rfloor$ 中 $0\sim m^2$ 是否出现过，以及需要计算前驱后继，可以做到 $O(m^2)$，所以总复杂度为 $O(\frac{p}{m}+m^2)$，把复杂度拆成 $\frac{p}{m\times2}+\frac{p}{m\times2}+m^2$，其几何平均数为 $\frac{p^{\frac23}}{2^{\frac23}}$，由均值不等式，我们可以得知总时间复杂度最低就是 $p^{\frac{2}{3}}$ 量级的，$m$ 取 $\frac{p^{\frac{1}{3}}}{2^{\frac13}}$ 时最优，实际取 $m=p^{\frac13}$ 即可，可根据实际情况调整。
自此，我们得到了 $O(p^{\frac{2}{3}})-O(1)$ 的优秀时间复杂度的在线逆元。实际表现拉完了，常数巨大，不过还是比快速幂和扩欧快一点，一般的数论题也不卡带 $\log$ 的逆元，这种算法适用于某些卡常题，正赛并不建议写。
同时，给大家一道例题测试一下自己的快速逆元吧！~~其实是测试卡常技巧。~~
:::success[[例题](/problem/P5431) && Code]
```cpp
#include<bits/stdc++.h>
#include<sys/mman.h>
#include<sys/stat.h>
#include<fcntl.h>
#include<unistd.h>
#define int long long
using namespace std;
const int maxn=1e6+7;
int inv[maxn],n,m,mm,pm,p,pre[maxn],nxt[maxn];
pair<int,int>f[maxn];
typedef unsigned long long ull;
ull M;
inline __attribute__((always_inline)) ull Mod(ull a){
ull r=a-(ull)((__int128)M*a>>64)*p;
return r>=p?r-p:r;
}
inline __attribute__((always_inline)) ull Div(ull a){
ull q=ull((__int128)M*a>>64);
ull r=a-q*p;
while(r<p)return q;
return q+1;
}
inline void init(){
m=pow(p,1.0/3);
mm=m*m;pm=p/m;
//mm 是 Farey 序列的项数，用于进行和 Farey 序列有关的预处理。
//pm 是 预处理逆元的项数，用于进行和逆元有关的计算。
//由于一些精度问题，mm 和 pm 不能一起算。
inv[1]=1;
f[0]={0,1};
f[mm]={1,1};
//f 数组用来存 Farey 序列的各数。
pre[mm]=nxt[mm]=mm;
//pre 和 nxt 数组用于存前驱后继。
for(int i=2;i<=pm;i++)inv[i]=Mod(inv[p%i]*(p-p/i));//线性求出 1 至 pm 的逆元。
for(int i=2;i<=m;i++){
for(int j=1;j<i;j++){//找 Farey 序列的每个数
int pos=j*mm/i;
if(!pre[pos]){//如果 Farey 序列里没有 i/j 就加入此数。
pre[pos]=nxt[pos]=pos;
f[pos]={j,i};
}
}
}
//直接在 Farey 序列中出现过的数都被录入了，接下来求每个可能值的前驱后继。
for(int i=1;i<=mm;i++)if(!pre[i])pre[i]=pre[i-1];
for(int i=mm;i;i--)if(!nxt[i])nxt[i]=nxt[i+1];
}
inline __attribute__((always_inline)) int get_val(int a,int id){
pair<int,int>d=f[id];
int u=a*d.second-p*d.first;
int v=u;u=abs(u);//我们需要验算左式的值和绝对值用于计算逆元。
while(u>pm)return -1;//如果绝对值大于阈值则不符合条件。
while(v<0)return p-Mod(d.second*inv[u]);//若左式的值为负则需要用 p 减一下以取负数。
return Mod(d.second*inv[u]);
}
int k,a,t;
uint8_t* str;
inline int read(){
int k=0;uint8_t c=(*str++);
while(!isdigit(c))c=(*str++);
while(isdigit(c))k=k*10+(c&15),c=(*str++);
return k;
}
signed main(){
struct stat ss;
fstat(0,&ss);
str=(uint8_t*)mmap(NULL,ss.st_size,1,2,0,0);
t=read();p=read();k=read();
M=((__int128)1<<64)/p;init();int ans=0,c=1;
while(t--){
c=Mod(c*k);
a=read();
int anss,pos;
if(a<=pm)ans=Mod(ans+c*inv[a]);
else{
pos=Div(a*mm);
anss=get_val(a,pre[pos]);
if(anss==-1)anss=get_val(a,nxt[pos]);
ans=Mod(ans+c*anss);
}
}
cout<<ans;
return 0;
}
```
最慢点 362ms，由于卡了常，所以码风会比较怪异，其实注释还挺清楚吧。
:::
由于预处理带来的时空开销，这个方法一般适用于模数不太大的情况，一般的常见模数（如 $998244353$ 和 $10^9+7$）都可以用这种方式求逆元。
本文之所以没什么例题，是因为题目都是推完式子再去计算的，大部分题的难点不在求逆元上。如果线练练手可以去找一些推式子题做。