多项式算法的常数问题

## 引言 最近几天自学了多项式全家桶。  众所周知，多项式算法的常数都比较巨大。比如多项式多点求值，$n=64000$，时间复杂度$O(n\log^2 n)$，时限竟然开到了3000ms，导致直接暴力卡卡常都能水过去。（卡常卡的好累） WC2019讲课人JOHNKRAM曾经说过，多项式的题目，只需要让暴力跑不过去就够了。 这一篇文章就主要讲述一下多项式算法的常数问题。 因为多项式全家桶已（wo）经（shi）在（zai）日（lan）报（de）里（xie）了，这篇文章几乎不会提及多项式算法的推导过程，所以……不会多项式全家桶的同学请参考[这篇洛谷日报](https://www.luogu.org/blog/user7035/duo-xiang-shi-zong-jie) 先贴一下我的多项式全家桶[代码](https://www.luogu.org/paste/5h43zj92)： ### FFT、NTT、FWT 这三种变换可以说是最基础的多项式算法了。 题目链接：[FFT](https://www.luogu.org/problemnew/show/P3803)，[FWT](https://www.luogu.org/problemnew/show/P4717)。 ```cpp void FFT(complex*A,int type) { for(int i=0;i<limit;i++) if(i<r[i])swap(A[i],A[r[i]]); for(int mid=1;mid<limit;mid<<=1) { complex Wn(cos(Pi/mid),type*sin(Pi/mid)); for(int R=mid<<1,j=0;j<limit;j+=R) { complex w(1,0); for(int k=0;k<mid;k++,w=w*Wn) { complex x=A[j+k],y=w*A[j+mid+k]; A[j+k]=x+y; A[j+mid+k]=x-y; } } } if(type==-1) { for(int i=0;i<limit;i++) A[i]=A[i]/limit; } } ll omega[2][21][1<<20]; void pre()//预处理单位根，可以显著优化NTT常数 { for(int mid=1,i=0;i<=20;mid<<=1,i++) { ll w1=quick_pow(3,(MOD-1)/(mid<<1)); ll w2=quick_pow(3,MOD-1-(MOD-1)/(mid<<1)); omega[0][i][0]=omega[1][i][0]=1; for(int k=1;k<mid;k++) { omega[1][i][k]=omega[1][i][k-1]*w1%MOD; omega[0][i][k]=omega[0][i][k-1]*w2%MOD; } } } void NTT(ll*A,int type) { for(int i=0;i<limit;i++) if(i<r[i])swap(A[i],A[r[i]]); if(type==-1)type=0; for(int mid=1,i=0;mid<limit;mid<<=1,i++) { for(int R=mid<<1,j=0;j<limit;j+=R) { for(int k=0;k<mid;k++) { ll x=A[j+k]%MOD,y=omega[type][i][k]*A[j+mid+k]%MOD; A[j+k]=x+y; A[j+mid+k]=x-y; } } } for(int i=0;i<limit;i++) { A[i]%=MOD; if(A[i]<0)A[i]+=MOD; } if(type==0) { ll inv=quick_pow(limit,MOD-2); for(int i=0;i<limit;i++)A[i]=A[i]*inv%MOD; } } void FWT_or(ll*A,int type) { for(int mid=1;mid<limit;mid<<=1) for(int R=mid<<1,j=0;j<limit;j+=R) for(int k=0;k<mid;k++) if(type==1)A[j+mid+k]+=A[j+k]; else A[j+mid+k]-=A[j+k]; for(int i=0;i<limit;i++) { A[i]%=MOD; if(A[i]<0)A[i]+=MOD; } } void FWT_and(ll*A,int type) { for(int mid=1;mid<limit;mid<<=1) for(int R=mid<<1,j=0;j<limit;j+=R) for(int k=0;k<mid;k++) if(type==1)A[j+k]+=A[j+mid+k]; else A[j+k]-=A[j+mid+k]; for(int i=0;i<limit;i++) { A[i]%=MOD; if(A[i]<0)A[i]+=MOD; } } void FWT_xor(ll*A,int type) { for(int mid=1;mid<limit;mid<<=1) for(int R=mid<<1,j=0;j<limit;j+=R) for(int k=0;k<mid;k++) { ll x=A[j+k],y=A[j+mid+k]; A[j+k]=x+y; A[j+mid+k]=x-y; if(A[j+k]>MOD)A[j+k]-=MOD; if(A[j+mid+k]<0)A[j+mid+k]+=MOD; if(type==-1) { A[j+k]=A[j+k]*inv2%MOD; A[j+mid+k]=A[j+mid+k]*inv2%MOD; } } } ``` 为了方便，我们假设这三种变换的常数都是$1$，也就是$n\log n$。（当然这些算法的常数也是有的，而且还不小） ### 任意模数FFT（MTT） [题目链接](https://www.luogu.org/problemnew/show/P4245) 共轭优化后，其实本质上就是4次DFT啦：（目前这份代码无论是时间上还是空间上排名都十分靠前，请自动忽略指令集优化的FFT） ```cpp void Mul(ll*A,ll*B,ll*C,ll MOD,complex*a,complex*b) { complex da1,db1,dc1,dd1,da2,db2,dc2,dd2,Da1,Db1,Dc1,Dd1,Da2,Db2,Dc2,Dd2; for(int i=0;i<limit;i++)a[i]=complex(A[i]&M,A[i]>>15); for(int i=0;i<limit;i++)b[i]=complex(B[i]&M,B[i]>>15); FFT(a,1); FFT(b,1); for(int i=0;i<=(limit>>1);i++) { int j=(limit-i)&(limit-1); da1=(a[i]+a[j].conj())*complex(0.5,0); db1=(a[i]-a[j].conj())*complex(0,-0.5); dc1=(b[i]+b[j].conj())*complex(0.5,0); dd1=(b[i]-b[j].conj())*complex(0,-0.5); da2=(a[j]+a[i].conj())*complex(0.5,0); db2=(a[j]-a[i].conj())*complex(0,-0.5); dc2=(b[j]+b[i].conj())*complex(0.5,0); dd2=(b[j]-b[i].conj())*complex(0,-0.5); Da1=da1*dc1; Db1=da1*dd1; Dc1=db1*dc1; Dd1=db1*dd1; Da2=da2*dc2; Db2=da2*dd2; Dc2=db2*dc2; Dd2=db2*dd2; a[i]=Da2+Db2*complex(0,1.0); b[i]=Dc2+Dd2*complex(0,1.0); a[j]=Da1+Db1*complex(0,1.0); b[j]=Dc1+Dd1*complex(0,1.0); } FFT(a,1); FFT(b,1); for(int i=0;i<limit;i++) { ll da=(ll)(a[i].x/limit+0.5)%MOD; ll db=(ll)(a[i].y/limit+0.5)%MOD; ll dc=(ll)(b[i].x/limit+0.5)%MOD; ll dd=(ll)(b[i].y/limit+0.5)%MOD; C[i]=(da+((ll)(db+dc)<<15)+((ll)dd<<30))%MOD; if(C[i]<0)C[i]+=MOD; } } ``` $$T(n)=4n\log n$$ ### 牛顿迭代 牛顿迭代计算复杂度的公式是 $$T(n)=T\left(n\over 2\right)+O(n\log n)=O(n\log n)$$ 让我们认真化一下计算复杂度的式子（我们先假设后面那个$O(n\log n)$的常数为$1$，所有的$\log$都是以$2$为底）： $$T(n)=T\left(n\over 2\right)+n\log n$$ $$=n\log n+\dfrac n2 \log\left(\dfrac n2\right)+T\left(n\over 4\right)$$ $$=\cdots=\sum_{i=0}^{\infty}\dfrac{n}{2^i}\log\left(\dfrac{n}{2^i}\right)$$ $$=\sum_{i=0}^{\infty}\dfrac{n}{2^i}\left(\log n-\log 2^i\right)$$ $$=n\log n\left(\sum_{i=0}^{\infty}\dfrac{1}{2^i}\right)-n\sum_{i=0}^{\infty}\dfrac{i}{2^i}$$ $$=2n\log n-2n$$ $$\approx 2n\log n$$ 也就是说，每一次迭代会让后面那个$n\log n$常数乘$2$。(为了方便我们忽略掉一次项） 利用这个式子，我们就可以得到其它多项式算法的常数： ### 多项式求逆 [题目链接](https://www.luogu.org/problemnew/show/P4238) 多项式求逆的原理是牛顿迭代+多项式乘法。 ```cpp void Inv(ll*a,ll*b,int n) { if(n==1) { b[0]=quick_pow(a[0],MOD-2); return; } Inv(a,b,(n+1)>>1);//牛顿迭代 limit=1,l=0; while(limit<(n<<1))limit<<=1,l++; for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); for(int i=0;i<limit;i++)A[i]=i<n?a[i]:0; for(int i=(n+1)>>1;i<limit;i++)b[i]=0; NTT(A,1);//nlogn NTT(b,1);//nlogn for(int i=0;i<limit;i++) b[i]=b[i]*(2+MOD-A[i]*b[i]%MOD)%MOD; NTT(b,-1);//nlogn for(int i=n;i<limit;i++)b[i]=0; } ``` $$T(n)=2T\left(\dfrac{n}{2}\right)+3n\log n=6n\log n$$ 我们可以看到，仅仅一个多项式求逆，常数就是NTT的6倍，多项式算法常数之大也就可想而知了。 ### 多项式开方 [题目链接](https://www.luogu.org/problemnew/show/P5205) 在多项式求逆的基础上，再进行一次牛顿迭代，想想就害怕。 ```cpp void Sqrt(ll*a,ll*b,int n) { if(n==1) { b[0]=1; return; } Sqrt(a,b,(n+1)>>1);//牛顿迭代 Inv(b,inv,n);//6nlogn limit=1,l=0; while(limit<(n<<1))limit<<=1,l++; for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); for(int i=0;i<limit;i++) A[i]=i<n?a[i]:0; NTT(inv,1);//nlogn NTT(A,1);//nlogn for(int i=0;i<limit;i++)inv[i]=inv[i]*A[i]%MOD; NTT(inv,-1);//nlogn for(int i=0;i<limit;i++)b[i]=(b[i]+inv[i])*inv2%MOD; for(int i=n;i<limit;i++)b[i]=0; } ``` $$T(n)=2T\left(\dfrac{n}{2}\right)+6n\log n+3n\log n=18n\log n$$ 常数已经上天了。 ### 多项式求ln： [题目链接](https://www.luogu.org/problemnew/show/P4725) 这个算法不需要牛顿迭代，但是要调用一次多项式求逆和三次NTT： ```cpp void Ln(ll*a,ll*b,int n) { Inv(a,b,n);//6nlogn limit=1,l=0; while(limit<(n<<1))limit<<=1,l++; for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); for(int i=0;i<n-1;i++)a2[i]=a[i+1]*(i+1)%MOD; for(int i=n-1;i<limit;i++)a2[i]=0; NTT(a2,1);//nlogn NTT(b,1);//nlogn for(int i=0;i<limit;i++) b[i]=b[i]*a2[i]%MOD; NTT(b,-1);//nlogn for(int i=n-1;i>0;i--) b[i]=b[i-1]*ni[i]%MOD; for(int i=n;i<limit;i++)b[i]=0; b[0]=0; } ``` $$T(n)=9n\log n$$ ### 多项式exp： [题目链接](https://www.luogu.org/problemnew/show/P4726) 先求多项式ln，然后牛顿迭代： ```cpp void Exp(ll*a,ll*b,int n) { if(n==1) { b[0]=1; return; } Exp(a,b,(n+1)>>1);//牛顿迭代 Ln(b,lnb,n);//9nlogn limit=1,l=0; while(limit<(n<<1))limit<<=1,l++; for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); for(int i=0;i<n;i++)lnb[i]=a[i]>=lnb[i]?a[i]-lnb[i]:a[i]-lnb[i]+MOD; for(int i=n;i<limit;i++)lnb[i]=b[i]=0; lnb[0]++; NTT(lnb,1);//nlogn NTT(b,1);//nlogn for(int i=0;i<limit;i++)b[i]=b[i]*lnb[i]%MOD; NTT(b,-1);//nlogn for(int i=n;i<limit;i++)b[i]=0; } ``` $$T(n)=2T\left(\dfrac{n}{2}\right)+9n\log n+3n\log n=24n\log n$$ （哪些多项式exp跑的比我的NTT还快的，你们是魔鬼吗） ### 多项式三角函数 [题目链接](https://www.luogu.org/problemnew/show/P5264) 本质上就是两遍多项式exp啦： ```cpp void Sin(ll*a,ll*b,int n) { for(int i=0;i<n;i++)a[i]=a[i]*img%MOD; Exp(a,tmp3,n);//24nlogn for(int i=0;i<n;i++)a[i]=MOD-a[i]; Exp(a,b,n);//24nlogn ll inv=inv2*invi%MOD; for(int i=0;i<n;i++) b[i]=(tmp3[i]-b[i]+MOD)*inv%MOD; } void Cos(ll*a,ll*b,int n) { for(int i=0;i<n;i++)a[i]=a[i]*img%MOD; Exp(a,tmp3,n);//24nlogn for(int i=0;i<n;i++)a[i]=MOD-a[i]; Exp(a,b,n);//24nlogn for(int i=0;i<n;i++) b[i]=(tmp3[i]+b[i])*inv2%MOD; } ``` $$T(n)=48n\log n$$ ### 多项式求arcsin [题目链接](https://www.luogu.org/problemnew/show/P5265) 查导数表可知： $$\arcsin'(x)=\dfrac{1}{\sqrt{1-x^2}}$$ 多项式平方+多项式开方+多项式求逆+多项式乘法： ```cpp void Asin(ll*f,ll*g,int n) { for(int i=1;i<n;i++) tmp3[i-1]=f[i]*i%MOD; limit=1,l=0; while(limit<(n<<1))limit<<=1,l++; for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); NTT(f,1);//nlogn for(int i=0;i<limit;i++)f[i]=f[i]*f[i]%MOD; NTT(f,-1);//nlogn for(int i=n;i<limit;i++)f[i]=0; for(int i=0;i<n;i++)f[i]=MOD-f[i]; f[0]++; Sqrt(f,g,n);//18nlogn Inv(g,f,n);//6nlogn limit=1,l=0; while(limit<(n<<1))limit<<=1,l++; for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); NTT(f,1);//nlogn NTT(tmp3,1);//nlogn for(int i=0;i<limit;i++)f[i]=f[i]*tmp3[i]%MOD; NTT(f,-1);//nlogn for(int i=n-1;i>=0;i--)g[i]=f[i-1]*ni[i]%MOD; g[0]=0; } ``` $$T(n)=29n\log n$$ ### 多项式求arctan 一样查导数表： $$\arctan'(x)=\dfrac{1}{1+x^2}$$ 多项式平方+多项式求逆+多项式乘法： ```cpp void Atan(ll*f,ll*g,int n) { for(int i=1;i<n;i++) tmp2[i-1]=f[i]*i%MOD; limit=1,l=0; while(limit<(n<<1))limit<<=1,l++; for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); NTT(f,1);//nlogn for(int i=0;i<limit;i++)f[i]=f[i]*f[i]%MOD; NTT(f,-1);//nlogn for(int i=n;i<limit;i++)f[i]=0; f[0]++; Inv(f,g,n);//6nlogn limit=1,l=0; while(limit<(n<<1))limit<<=1,l++; for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); NTT(g,1);//nlogn NTT(tmp2,1);//nlogn for(int i=0;i<limit;i++)g[i]=g[i]*tmp2[i]%MOD; NTT(g,-1);//nlogn for(int i=n-1;i>=0;i--)g[i]=g[i-1]*ni[i]%MOD; g[0]=0; } ``` 常数稍小一点， $$T(n)=11n\log n$$ ### 多项式带余除法、多项式取模 [题目链接](https://www.luogu.org/problemnew/show/P4512) 一遍求逆+6遍NTT： ```cpp void Div(ll*F,ll*G,ll*Q,ll*R,int n,int m) { for(int i=0;i<(m>>1);i++)swap(G[i],G[m-i-1]); Inv(G,Q,n-m+1);//6nlogn limit=1,l=0; while(limit<(n<<1))limit<<=1,l++; for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); for(int i=0;i<n-m+1;i++)tmp1[i]=F[n-i-1]; for(int i=n-m+1;i<limit;i++)tmp1[i]=0; NTT(Q,1);//nlogn NTT(tmp1,1);//nlogn for(int i=0;i<limit;i++)Q[i]=Q[i]*tmp1[i]%MOD; NTT(Q,-1);//nlogn for(int i=0;i<(m>>1);i++)swap(G[i],G[m-i-1]); for(int i=0;i<((n-m+1)>>1);i++)swap(Q[i],Q[n-m-i]); for(int i=n-m+1;i<limit;i++)Q[i]=0; for(int i=0;i<limit;i++)tmp1[i]=Q[i],R[i]=i<m?G[i]:0; NTT(tmp1,1);//nlogn NTT(R,1);//nlogn for(int i=0;i<limit;i++)R[i]=tmp1[i]*R[i]%MOD; NTT(R,-1);//nlogn for(int i=0;i<m-1;i++)R[i]=(F[i]>=R[i]?F[i]-R[i]:F[i]-R[i]+MOD); for(int i=m-1;i<limit;i++)R[i]=0; } void Mod(ll*F,ll*G,ll*R,int n,int m) { if(n<m) { for(int i=0;i<n;i++)R[i]=F[i]; return; } Div(F,G,tmp2,R,n,m); } ``` $$T(n)=12n\log n$$ ### 主定理的常数 再往下，时间复杂度终于不是$O(n\log n)$了。 由主定理可知， $$T(n)=2T\left(\dfrac n2\right)+O(n\log n)\implies T(n)=O(n\log^2 n)$$ 那么主定理的常数又有多少呢？ 一个$T(n)$会变成两个$T\left(\dfrac n2\right)$，所以我们可以将复杂度画成一个树形结构：  容易发现最终的复杂度就是最后一列所有的数加起来。 $$\sum_{i=0}^{\log n}n\log \dfrac{n}{2^i}$$ $$=n\sum_{i=0}^{\log n}(\log n-\log 2^i)$$ $$=n\log^2 n-n\sum_{i=0}^{\log n}i$$ $$=n\log^2 n-n\dfrac{\log n(\log n+1)}{2}$$ $$\approx\dfrac{1}{2}n\log^2 n$$ 也就是说，一次主定理会让复杂度多一个$\log$，但是常数减半。注意这里的常数如果转化到一个$\log$的话大概还要乘$20$左右。 ### 多项式多点求值 [题目链接](https://www.luogu.org/problemnew/show/P5050) 多项式多点求值的思路是什么呢？ $$f(x_i)=f(x)\bmod{\:(x-x_i)}$$ 我们需要先用分治求出$\prod_{i=l}^r(x-x_i)$，然后再分治取模。 这样的时间复杂度分为两部分： ```cpp void solve(int o,int L,int R,ll*x) { if(L==R) { _p[o]=__p+tot; tmp3[o]=2; tot+=2; _p[o][0]=MOD-x[L]; _p[o][1]=1; return; } int mid=(L+R)>>1; solve(o<<1,L,mid,x); solve(o<<1|1,mid+1,R,x);//2T(n/2) limit=1,l=0; while(limit<=(R-L+1))limit<<=1,l++; for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); for(int i=0;i<tmp3[o<<1];i++)tmp1[i]=_p[o<<1][i]; for(int i=tmp3[o<<1];i<limit;i++)tmp1[i]=0; for(int i=0;i<tmp3[o<<1|1];i++)tmp2[i]=_p[o<<1|1][i]; for(int i=tmp3[o<<1|1];i<limit;i++)tmp2[i]=0; NTT(tmp1,1);//nlogn NTT(tmp2,1);//nlogn for(int i=0;i<limit;i++)tmp1[i]=tmp1[i]*tmp2[i]%MOD; NTT(tmp1,-1);//nlogn while(limit>1&&tmp1[limit-1]==0)limit--; _p[o]=__p+tot; tmp3[o]=limit; tot+=limit; for(int i=0;i<limit;i++)_p[o][i]=tmp1[i]; } void query(int o,int dep,int l,int r,ll*ans) { if(l==r) { ans[l]=_f[dep][0]; return; } int mid=(l+r)>>1; Mod(_f[dep],_p[o<<1],_f[dep+1],tmp3[o],tmp3[o<<1]);//12nlogn query(o<<1,dep+1,l,mid,ans);//T(n/2) Mod(_f[dep],_p[o<<1|1],_f[dep+1],tmp3[o],tmp3[o<<1|1]);//12nlogn query(o<<1|1,dep+1,mid+1,r,ans);//T(n/2) } void Eval(ll*a,ll*x,int n,int m,ll*res) { solve(1,1,m,x); Mod(a,_p[1],_f[1],n,tmp3[1]); query(1,1,1,m,res); } ``` $$T(n)=T_1(n)+T_2(n)=\left(2T_1\left(\dfrac n2\right)+3n\log n\right)+\left(2T_2\left(\dfrac n2\right)+24n\log n\right)$$ $$T(n)=1.5n\log^2 n+12n\log^2 n=13.5n\log^2 n$$ $n=64000$，时限3000ms，就是这么来的，注意这里还忽略了$n\log n$的项。 ### 多项式快速插值 （啥快速啊，没准还不如暴力快） [题目链接](https://www.luogu.org/problemnew/show/P5158) 多项式快速插值的原理是拉格朗日插值公式，也就是 $$A(x)=\sum_{i=1}^ny_i\dfrac{\prod_{j\neq i}(x-x_j)}{\prod_{j\neq i}(x_i-x_j)}$$ 注意到$\prod_{j\neq i}(x-x_j)=\pi'(x_i)$，其中$\pi(x)=\prod_{i=1}^n(x-x_i)$。 所以我们先用分治乘法求出$\pi(x)$，求导后多点求值，最后再用一个分治乘法求出分母部分。 这样做的常数问题嘛……还是先上代码： ```cpp void ask(int o,int dep,int L,int R,ll*y) { if(L==R) { _f[dep][0]=y[L]; return; } int mid=(L+R)>>1; ask(o<<1,dep+1,L,mid,y);//T(n/2) limit=1,l=0; while(limit<R-L+1)limit<<=1,l++; for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); for(int i=0;i<mid-L+1;i++)tmp1[i]=_f[dep+1][i]; for(int i=mid-L+1;i<limit;i++)tmp1[i]=0; for(int i=0;i<tmp3[o<<1|1];i++)tmp2[i]=_p[o<<1|1][i]; for(int i=tmp3[o<<1|1];i<limit;i++)tmp2[i]=0; NTT(tmp1,1);//nlogn NTT(tmp2,1);//nlogn for(int i=0;i<limit;i++)_f[dep][i]=tmp1[i]*tmp2[i]; ask(o<<1|1,dep+1,mid+1,R,y);//T(n/2) limit=1,l=0; while(limit<R-L+1)limit<<=1,l++; for(int i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1)); for(int i=0;i<tmp3[o<<1];i++)tmp1[i]=_p[o<<1][i]; for(int i=tmp3[o<<1];i<limit;i++)tmp1[i]=0; for(int i=0;i<R-mid;i++)tmp2[i]=_f[dep+1][i]; for(int i=R-mid;i<limit;i++)tmp2[i]=0; NTT(tmp1,1);//nlogn NTT(tmp2,1);//nlogn for(int i=0;i<limit;i++)_f[dep][i]=(_f[dep][i]+tmp1[i]*tmp2[i])%MOD; NTT(_f[dep],-1);//nlogn } void Inter(ll*x,ll*y,int n,ll*res) { solve(1,1,n,x);//求出\pi(x) for(int i=1;i<tmp3[1];i++) res[i-1]=_p[1][i]*i%MOD; for(int i=tmp3[1]-1;i<tmp3[1]<<2;i++)res[i]=0; Mod(res,_p[1],_f[1],tmp3[1]-1,tmp3[1]); query(1,1,1,n,res);//带入多点求值 for(int i=1;i<=n;i++)res[i]=y[i]*quick_pow(res[i],MOD-2)%MOD; ask(1,1,1,n,res);//最后再分治乘一波 for(int i=0;i<n;i++)res[i]=_f[1][i]%MOD; } ``` $$T(n)=T_1(n)+T_2(n)=13.5n\log^2 n+\left(T_2\left(\dfrac{n}{2}\right)+5n\log n\right)=16n\log^2 n$$ ## 总结 总之就是常数大。 由于忽略了所有的低阶项，所以这些算法的常数会比理论计算要大一些。 理论要有实践来支撑。在Linux虚拟机里面运行的结果如下（没有开O2，三次测量四舍五入取平均）： | 算法 | 理论复杂度 | $n=100000$ | 与NTT的比值 | $n=1000000$ | 与NTT的比值 | | :----------: | :----------: | :----------: | :----------: | :----------: | :----------: | | NTT | $n\log n$ | 50ms | 1 | 500ms | 1 | | FFT | $n\log n$ | 124ms | 2.48 | 900ms | 1.8 | | FWT(or) | $n\log n$ | 40ms | 0.8 | 170ms | 0.34 | | FWT(and) | $n\log n$ | 50ms | 1 | 175ms | 0.35 | | FWT(xor) | $n\log n$ | 70ms | 1.4 | 330ms | 0.66 | | MTT | $4n\log n$ | 870ms | 17.4 | 7620ms | 15.24 | | 多项式求逆 | $6n\log n$ | 528ms | 10.56 | 4680ms | 9.36 | | 多项式开方 | $18n\log n$ | 1404ms | 28.08 | 13300ms | 26.6 | | 多项式$\ln$ | $9n\log n$ | 770ms | 15.4 | 7150ms | 14.3 | | 多项式$\exp$ | $24n\log n$ | 1860ms | 37.2 | 17800ms | 35.6 | | 多项式$\sin$ | $48n\log n$ | 3680ms | 73.6 | 35100ms | 70.2 | | 多项式$\cos$ | $48n\log n$ | 3690ms | 73.8 | 35200ms | 70.4 | | 多项式$\arcsin$ | $29n\log n$ | 2290ms | 45.8 | 21700ms | 43.4 | | 多项式$\arctan$ | $11n\log n$ | 940ms | 18.8 | 8650ms | 17.3 | | 多项式带余除法 | $12n\log n$ | 760ms | 15.2 | 7000ms | 14 | | 多项式多点求值 | $13.5n\log n$ | 15400ms | 308 | 175000ms | 350 | | 多项式快速插值 | $16n\log n$ | 17400ms | 348 | 196500ms | 393 | （你能想象3分钟一个点的**快速**插值吗） 由于虚拟机的运行速度可想而知，以上所有运行时间参考意义都不大，放到洛谷上大概要除以8（我的快速插值在洛谷上$100000$的数据最慢的点跑了2000ms） 所有运行时间小于1000ms的误差都相当大，而且我的NTT前面加了一个预处理的常数，所以前面FFT和FWT与NTT运行时间的比值会有较大的误差。 可以发现在实际运行中，所有算法与NTT的相对常数大概还要再乘$1.5$左右，其中任意模数NTT更是比理论复杂度慢了2倍（与FFT相比）。但是随着$n$的增大，时间复杂度相同的情况下，绝大多数比值在逐渐减小，这也符合渐进复杂度的规律。 ## 多项式常数大，慎用！