其实一直想做这个测试了,但是没时间(懒)。
某个下午没干劲做题,不想荒废于是搞了出来。
引子
这种求逆元的方法其实是半年前指挥安利给我的。
众所周知,逆元可以$O(n)$递推:
const int mod = 998244353;
int inv[maxn],n;
int main(){
for(register int i=1;i<=n;++i)inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
}
这只适用于模数较小或用的逆元个数有限的情况。
但我们可以把递推改成递归:
const int mod = 998244353;
int inv(int x){return x==1?1:1ll*(mod-mod/x)*inv(mod%x)%mod;}
这样可以求出任意数的逆元。如果会出现$0$要特判。
这种方法好像也没有什么专业的名称,我姑且称其为「递归式逆元」
优势非常明显,特别短,只有一行。
但它依然不适用于模数非质数的情况,而且复杂度未知。
接下来就是关于其复杂度的胡搞。
尝试证明
这个复杂度不好证明。
每次$x$会变为$mod\%x$,这并不能保证$x$会缩小为原来的多少。($exgcd$每次$a$至少缩小一半所以保证其复杂度是$\log$的)
在洛谷上得到了一篇知乎的证明:戳这儿。
大概意思是上界大约为$O(n^{\frac{1}{3}})$,下界为$O(\dfrac{\ln n}{\ln\ln n})$,期望为$O(\log mod)$。
这个复杂度听起来很爆炸。
好证明不如烂测试。直接挑选几个模数测试$1\sim mod-1$里最大递归次数:
const int mod = 998244353;
int cnt[maxn];
int inv(int x){if(x<maxn)return cnt[x];return inv(mod%x)+1;}
int main(){
int ans=0;
for(register int i=2;i<maxn;++i)ans=max(ans,cnt[i]=cnt[mod%i]+1);
for(register int i=maxn;i<mod;++i)ans=max(ans,inv(i));
printf("%d\n",ans);
}
结果:
模数 | 最大递归次数 |
---|---|
$998244353$ | $44$ |
$10^9+7$ | $48$ |
$10^9+9$ | $46$ |
$10^6+7$ | $25$ |
$104857601$ | $38$ |
某$8$位大质数 | $36$ |
这至少说明了在一些常用模数下,复杂度还是可以接受的。
可以看出最大递归次数与模数是近似的正相关,但是也有$10^9+7$和$10^9+9$的反例。
比较
其他求逆元的方法无非就是基于费马小定理的快速幂和$exgcd$了。
论适用性的话还是$exgcd$,能处理模数非质数的情况。
以下只比较效率。
由于取模和除法效率远低于其他运算,只考虑这两种。同时还考虑了递归开栈和传入变量。
快速幂
inline int quickpow(int x,int y){
int ans=1;
while(y){
if(y&1)ans=1ll*ans*x%mod;
x=1ll*x*x%mod,y>>=1;
}
}
运行$\log mod$次,一般跑不满。每次涉及$2$次取模,没有递归。
exgcd
void exgcd(int a,int b,int &x,int &y){
if(!b)x=1,y=0;
else exgcd(b,a%b,y,x),y-=x*(a/b);
}
inline int inv(int a){
int x,y;
exgcd(a, mod, x, y);
return (x+mod)%mod;
}
运行$\log mod$次,一般跑不满。每次涉及$1$次取模、$1$次除法。
需要递归实现,每层递归传入$2$个$int$和$2$个$int\&$。
递归式
int inv(int x){return x==1?1:1ll*(mod-mod/x)*inv(mod%x)%mod;}
运行$?$次(参考上面数据)。每次涉及$2$次取模、$1$次除法。
需要递归实现,每层传入$1$个$int$。
小结
理论上快速幂的性能是最高的,$exgcd$常数不小,递归式复杂度和常数都很大。
但是前两种递归式都不好优化,递归式有一个非常棒的优化:预处理。
int inv[maxn]={0,1};
int INV(int x){
if(x<maxn)return inv[x];
return 1ll*(mod-mod/x)*INV(mod%x)%mod;
}
int main(){
for(register int i=2;i<maxn;++i)inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
}
也可以记忆化,但估计是负优化而且难写。
然而以上都是纸上谈兵,只是理论瞎分析了一波常数。实测才是硬道理。
实测
以下测试都只用最常见的$998244353$和$10^9+7$,因为clock()
是厌氧生物所以都没开$O2$。
当初指挥给我安利递归式的时候,进行了随机数据下的测试,结果是递归式跑的飞起。
那么把$1\sim mod-1$全跑一遍呢?
对快速幂、$exgcd$、递归式、预处理递归式(预处理至$10^7$,占用空间$38MB$,测试时算上预处理时间)进行测试。
测试代码:
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <time.h>
#define maxn 10000005
const int mod = 998244353;
using namespace std;
namespace inv1{
inline int quickpow(int x,int y){
int ans=1;
while(y){
if(y&1)ans=1ll*ans*x%mod;
x=1ll*x*x%mod,y>>=1;
}
return ans;
}
}
namespace inv2{
void exgcd(int a,int b,int &x,int &y){
if(!b)x=1,y=0;
else exgcd(b,a%b,y,x),y-=x*(a/b);
}
inline int inv(int a){
int x,y;
exgcd(a,mod,x,y);
return (x+mod)%mod;
}
}
namespace inv3{
int inv(int x){return x==1?1:1ll*(mod-mod/x)*inv(mod%x)%mod;}
}
namespace inv4{
int inv[maxn];
void sieve(){
inv[1]=1;
for(register int i=2;i<maxn;++i)inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
}
int INV(int x){
if(x<maxn)return inv[x];
return x==1?1:1ll*(mod-mod/x)*INV(mod%x)%mod;
}
}
int main(){
printf("Speed test under mod %d.\n",mod);
int s=clock();
for(register int i=1;i<mod;++i)inv1::quickpow(i,mod-2);
s=clock()-s;
printf("quickpow : %d ms\n",s);
s=clock();
for(register int i=1;i<mod;++i)inv2::inv(i);
s=clock()-s;
printf("exgcd : %d ms\n",s);
s=clock();
for(register int i=1;i<mod;++i)inv3::inv(i);
s=clock()-s;
printf("recursive : %d ms\n",s);
s=clock();
for(register int i=1;i<mod;++i)inv4::INV(i);
s=clock()-s;
printf("recursive with preprocess: %d ms\n",s);
}
不要在意我胡写的英语
结果:
$998244353$:
方法 | 时间(ms) |
---|---|
快速幂 | 162268 |
$exgcd$ | 260005 |
递归式 | 270803 |
预处理递归式 | 99148 |
$10^9+7$也没啥区别:
方法 | 时间(ms) |
---|---|
快速幂 | 175463 |
$exgcd$ | 270038 |
递归式 | 328440 |
预处理递归式 | 103875 |
预处理竟恐怖如斯,甚至一度突破$100s$。
缺点在于占用空间。可以少预处理一些以时间换空间。
再一个预处理递归式适用范围可能不是很广,并没有什么题会让你把$1\sim mod-1$的逆元全求一遍。拿来卡常的话预处理自身已经需要$200+ms$了。
其他方法里快速幂是最优秀的了,而$exgcd$和递归式常数确实大。
总结
逆元用什么方法还是看习惯吧。
通过这些测试,我收获了:
- 一篇水文
- 颓废的快乐
- 少量的干劲
- “没有荒废”的自我安慰
以上。