多项式算法的常数问题

引言

最近几天自学了多项式全家桶。

众所周知，多项式算法的常数都比较巨大。比如多项式多点求值，n=64000n=64000n=64000，时间复杂度O(nlog⁡2n)O(n\log^2 n)O(nlog2n)，时限竟然开到了3000ms，导致直接暴力卡卡常都能水过去。（卡常卡的好累）

WC2019讲课人JOHNKRAM曾经说过，多项式的题目，只需要让暴力跑不过去就够了。

这一篇文章就主要讲述一下多项式算法的常数问题。

因为多项式全家桶已（wo）经（shi）在（zai）日（lan）报（de）里（xie）了，这篇文章几乎不会提及多项式算法的推导过程，所以……不会多项式全家桶的同学请参考这篇洛谷日报

先贴一下我的多项式全家桶代码：

FFT、NTT、FWT

这三种变换可以说是最基础的多项式算法了。

题目链接：FFT，FWT。

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;}}
}

为了方便，我们假设这三种变换的常数都是111，也就是nlog⁡nn\log nnlogn。（当然这些算法的常数也是有的，而且还不小）

任意模数FFT（MTT）

题目链接

共轭优化后，其实本质上就是4次DFT啦：（目前这份代码无论是时间上还是空间上排名都十分靠前，请自动忽略指令集优化的FFT）

void Mul(ll*A,ll*B,ll*C,ll MOD,complex*a,complex*b)//4nlogn
{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)=4nlog⁡nT(n)=4n\log nT(n)=4nlogn

牛顿迭代

牛顿迭代计算复杂度的公式是

T(n)=T(n2)+O(nlog⁡n)=O(nlog⁡n)T(n)=T\left(n\over 2\right)+O(n\log n)=O(n\log n)T(n)=T(2n)+O(nlogn)=O(nlogn)

让我们认真化一下计算复杂度的式子（我们先假设后面那个O(nlog⁡n)O(n\log n)O(nlogn)的常数为111，所有的log⁡\loglog都是以222为底）：

T(n)=T(n2)+nlog⁡nT(n)=T\left(n\over 2\right)+n\log nT(n)=T(2n)+nlogn

=nlog⁡n+n2log⁡(n2)+T(n4)=n\log n+\dfrac n2 \log\left(\dfrac n2\right)+T\left(n\over 4\right)=nlogn+2nlog(2n)+T(4n)

=⋯=∑i=0∞n2ilog⁡(n2i)=\cdots=\sum_{i=0}^{\infty}\dfrac{n}{2^i}\log\left(\dfrac{n}{2^i}\right)=⋯=i=0∑∞2inlog(2in)

=∑i=0∞n2i(log⁡n−log⁡2i)=\sum_{i=0}^{\infty}\dfrac{n}{2^i}\left(\log n-\log 2^i\right)=i=0∑∞2in(logn−log2i)

=nlog⁡n(∑i=0∞12i)−n∑i=0∞i2i=n\log n\left(\sum_{i=0}^{\infty}\dfrac{1}{2^i}\right)-n\sum_{i=0}^{\infty}\dfrac{i}{2^i}=nlogn(i=0∑∞2i1)−ni=0∑∞2ii

=2nlog⁡n−2n=2n\log n-2n=2nlogn−2n

≈2nlog⁡n\approx 2n\log n≈2nlogn

也就是说，每一次迭代会让后面那个nlog⁡nn\log nnlogn常数乘222。(为了方便我们忽略掉一次项）

利用这个式子，我们就可以得到其它多项式算法的常数：

多项式求逆

题目链接

多项式求逆的原理是牛顿迭代+多项式乘法。

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);//nlognNTT(b,1);//nlognfor(int i=0;i<limit;i++)b[i]=b[i]*(2+MOD-A[i]*b[i]%MOD)%MOD;NTT(b,-1);//nlognfor(int i=n;i<limit;i++)b[i]=0;
}

T(n)=2T(n2)+3nlog⁡n=6nlog⁡nT(n)=2T\left(\dfrac{n}{2}\right)+3n\log n=6n\log nT(n)=2T(2n)+3nlogn=6nlogn

我们可以看到，仅仅一个多项式求逆，常数就是NTT的6倍，多项式算法常数之大也就可想而知了。

多项式开方

题目链接

在多项式求逆的基础上，再进行一次牛顿迭代，想想就害怕。

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);//6nlognlimit=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);//nlognNTT(A,1);//nlognfor(int i=0;i<limit;i++)inv[i]=inv[i]*A[i]%MOD;NTT(inv,-1);//nlognfor(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(n2)+6nlog⁡n+3nlog⁡n=18nlog⁡nT(n)=2T\left(\dfrac{n}{2}\right)+6n\log n+3n\log n=18n\log nT(n)=2T(2n)+6nlogn+3nlogn=18nlogn

常数已经上天了。

多项式求ln：

题目链接

这个算法不需要牛顿迭代，但是要调用一次多项式求逆和三次NTT：

void Ln(ll*a,ll*b,int n)
{Inv(a,b,n);//6nlognlimit=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);//nlognNTT(b,1);//nlognfor(int i=0;i<limit;i++)b[i]=b[i]*a2[i]%MOD;NTT(b,-1);//nlognfor(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)=9nlog⁡nT(n)=9n\log nT(n)=9nlogn

多项式exp：

题目链接

先求多项式ln，然后牛顿迭代：

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);//9nlognlimit=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);//nlognNTT(b,1);//nlognfor(int i=0;i<limit;i++)b[i]=b[i]*lnb[i]%MOD;NTT(b,-1);//nlognfor(int i=n;i<limit;i++)b[i]=0;
}

T(n)=2T(n2)+9nlog⁡n+3nlog⁡n=24nlog⁡nT(n)=2T\left(\dfrac{n}{2}\right)+9n\log n+3n\log n=24n\log nT(n)=2T(2n)+9nlogn+3nlogn=24nlogn

（哪些多项式exp跑的比我的NTT还快的，你们是魔鬼吗）

多项式三角函数

题目链接

本质上就是两遍多项式exp啦：

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);//24nlognfor(int i=0;i<n;i++)a[i]=MOD-a[i];Exp(a,b,n);//24nlognll 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);//24nlognfor(int i=0;i<n;i++)a[i]=MOD-a[i];Exp(a,b,n);//24nlognfor(int i=0;i<n;i++)b[i]=(tmp3[i]+b[i])*inv2%MOD;
}

T(n)=48nlog⁡nT(n)=48n\log nT(n)=48nlogn

多项式求arcsin

题目链接

查导数表可知：

arcsin⁡′(x)=11−x2\arcsin'(x)=\dfrac{1}{\sqrt{1-x^2}}arcsin′(x)=1−x21

多项式平方+多项式开方+多项式求逆+多项式乘法：

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);//nlognfor(int i=0;i<limit;i++)f[i]=f[i]*f[i]%MOD;NTT(f,-1);//nlognfor(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);//18nlognInv(g,f,n);//6nlognlimit=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);//nlognNTT(tmp3,1);//nlognfor(int i=0;i<limit;i++)f[i]=f[i]*tmp3[i]%MOD;NTT(f,-1);//nlognfor(int i=n-1;i>=0;i--)g[i]=f[i-1]*ni[i]%MOD;g[0]=0;
}

T(n)=29nlog⁡nT(n)=29n\log nT(n)=29nlogn

多项式求arctan

一样查导数表：

arctan⁡′(x)=11+x2\arctan'(x)=\dfrac{1}{1+x^2}arctan′(x)=1+x21

多项式平方+多项式求逆+多项式乘法：

void Atan(ll*f,ll*g,int n)//11nlogn//tmp1,tmp2
{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);//nlognfor(int i=0;i<limit;i++)f[i]=f[i]*f[i]%MOD;NTT(f,-1);//nlognfor(int i=n;i<limit;i++)f[i]=0;f[0]++;Inv(f,g,n);//6nlognlimit=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);//nlognNTT(tmp2,1);//nlognfor(int i=0;i<limit;i++)g[i]=g[i]*tmp2[i]%MOD;NTT(g,-1);//nlognfor(int i=n-1;i>=0;i--)g[i]=g[i-1]*ni[i]%MOD;g[0]=0;
}

常数稍小一点，

T(n)=11nlog⁡nT(n)=11n\log nT(n)=11nlogn

多项式带余除法、多项式取模

题目链接

一遍求逆+6遍NTT：

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);//6nlognlimit=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);//nlognNTT(tmp1,1);//nlognfor(int i=0;i<limit;i++)Q[i]=Q[i]*tmp1[i]%MOD;NTT(Q,-1);//nlognfor(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);//nlognNTT(R,1);//nlognfor(int i=0;i<limit;i++)R[i]=tmp1[i]*R[i]%MOD;NTT(R,-1);//nlognfor(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)=12nlog⁡nT(n)=12n\log nT(n)=12nlogn

主定理的常数

再往下，时间复杂度终于不是O(nlog⁡n)O(n\log n)O(nlogn)了。

由主定理可知，

T(n)=2T(n2)+O(nlog⁡n) ⟹ T(n)=O(nlog⁡2n)T(n)=2T\left(\dfrac n2\right)+O(n\log n)\implies T(n)=O(n\log^2 n)T(n)=2T(2n)+O(nlogn)⟹T(n)=O(nlog2n)

那么主定理的常数又有多少呢？

一个T(n)T(n)T(n)会变成两个T(n2)T\left(\dfrac n2\right)T(2n)，所以我们可以将复杂度画成一个树形结构：

容易发现最终的复杂度就是最后一列所有的数加起来。

∑i=0log⁡nnlog⁡n2i\sum_{i=0}^{\log n}n\log \dfrac{n}{2^i}i=0∑lognnlog2in

=n∑i=0log⁡n(log⁡n−log⁡2i)=n\sum_{i=0}^{\log n}(\log n-\log 2^i)=ni=0∑logn(logn−log2i)

=nlog⁡2n−n∑i=0log⁡ni=n\log^2 n-n\sum_{i=0}^{\log n}i=nlog2n−ni=0∑logni

=nlog⁡2n−nlog⁡n(log⁡n+1)2=n\log^2 n-n\dfrac{\log n(\log n+1)}{2}=nlog2n−n2logn(logn+1)

≈12nlog⁡2n\approx\dfrac{1}{2}n\log^2 n≈21nlog2n

也就是说，一次主定理会让复杂度多一个log⁡\loglog，但是常数减半。注意这里的常数如果转化到一个log⁡\loglog的话大概还要乘202020左右。

多项式多点求值

题目链接

多项式多点求值的思路是什么呢？

f(xi)=f(x)&VeryThinSpace;mod&VeryThinSpace; (x−xi)f(x_i)=f(x)\bmod{\:(x-x_i)}f(xi)=f(x)mod(x−xi)

我们需要先用分治求出∏i=lr(x−xi)\prod_{i=l}^r(x-x_i)∏i=lr(x−xi)，然后再分治取模。

这样的时间复杂度分为两部分：

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);//nlognNTT(tmp2,1);//nlognfor(int i=0;i<limit;i++)tmp1[i]=tmp1[i]*tmp2[i]%MOD;NTT(tmp1,-1);//nlognwhile(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]);//12nlognquery(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]);//12nlognquery(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)=T1(n)+T2(n)=(2T1(n2)+3nlog⁡n)+(2T2(n2)+24nlog⁡n)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)=T1(n)+T2(n)=(2T1(2n)+3nlogn)+(2T2(2n)+24nlogn)

T(n)=1.5nlog⁡2n+12nlog⁡2n=13.5nlog⁡2nT(n)=1.5n\log^2 n+12n\log^2 n=13.5n\log^2 nT(n)=1.5nlog2n+12nlog2n=13.5nlog2n

n=64000n=64000n=64000，时限3000ms，就是这么来的，注意这里还忽略了nlog⁡nn\log nnlogn的项。

多项式快速插值

（啥快速啊，没准还不如暴力快）

题目链接

多项式快速插值的原理是拉格朗日插值公式，也就是

A(x)=∑i=1nyi∏j≠i(x−xj)∏j≠i(xi−xj)A(x)=\sum_{i=1}^ny_i\dfrac{\prod_{j\neq i}(x-x_j)}{\prod_{j\neq i}(x_i-x_j)}A(x)=i=1∑nyi∏j̸=i(xi−xj)∏j̸=i(x−xj)

注意到∏j≠i(x−xj)=π′(xi)\prod_{j\neq i}(x-x_j)=\pi'(x_i)∏j̸=i(x−xj)=π′(xi)，其中π(x)=∏i=1n(x−xi)\pi(x)=\prod_{i=1}^n(x-x_i)π(x)=∏i=1n(x−xi)。

所以我们先用分治乘法求出π(x)\pi(x)π(x)，求导后多点求值，最后再用一个分治乘法求出分母部分。

这样做的常数问题嘛……还是先上代码：

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);//nlognNTT(tmp2,1);//nlognfor(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);//nlognNTT(tmp2,1);//nlognfor(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)=T1(n)+T2(n)=13.5nlog⁡2n+(T2(n2)+5nlog⁡n)=16nlog⁡2nT(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 nT(n)=T1(n)+T2(n)=13.5nlog2n+(T2(2n)+5nlogn)=16nlog2n

总结

总之就是常数大。

由于忽略了所有的低阶项，所以这些算法的常数会比理论计算要大一些。

理论要有实践来支撑。在Linux虚拟机里面运行的结果如下（没有开O2，三次测量四舍五入取平均）：

算法	理论复杂度	n=100000n=100000n=100000	与NTT的比值	n=1000000n=1000000n=1000000	与NTT的比值
NTT	nlog⁡nn\log nnlogn	50ms	1	500ms	1
FFT	nlog⁡nn\log nnlogn	124ms	2.48	900ms	1.8
FWT(or)	nlog⁡nn\log nnlogn	40ms	0.8	170ms	0.34
FWT(and)	nlog⁡nn\log nnlogn	50ms	1	175ms	0.35
FWT(xor)	nlog⁡nn\log nnlogn	70ms	1.4	330ms	0.66
MTT	4nlog⁡n4n\log n4nlogn	870ms	17.4	7620ms	15.24
多项式求逆	6nlog⁡n6n\log n6nlogn	528ms	10.56	4680ms	9.36
多项式开方	18nlog⁡n18n\log n18nlogn	1404ms	28.08	13300ms	26.6
多项式ln⁡\lnln	9nlog⁡n9n\log n9nlogn	770ms	15.4	7150ms	14.3
多项式exp⁡\expexp	24nlog⁡n24n\log n24nlogn	1860ms	37.2	17800ms	35.6
多项式sin⁡\sinsin	48nlog⁡n48n\log n48nlogn	3680ms	73.6	35100ms	70.2
多项式cos⁡\coscos	48nlog⁡n48n\log n48nlogn	3690ms	73.8	35200ms	70.4
多项式arcsin⁡\arcsinarcsin	29nlog⁡n29n\log n29nlogn	2290ms	45.8	21700ms	43.4
多项式arctan⁡\arctanarctan	11nlog⁡n11n\log n11nlogn	940ms	18.8	8650ms	17.3
多项式带余除法	12nlog⁡n12n\log n12nlogn	760ms	15.2	7000ms	14
多项式多点求值	13.5nlog⁡n13.5n\log n13.5nlogn	15400ms	308	175000ms	350
多项式快速插值	16nlog⁡n16n\log n16nlogn	17400ms	348	196500ms	393

（你能想象3分钟一个点的快速插值吗）

由于虚拟机的运行速度可想而知，以上所有运行时间参考意义都不大，放到洛谷上大概要除以8（我的快速插值在洛谷上100000100000100000的数据最慢的点跑了2000ms）

所有运行时间小于1000ms的误差都相当大，而且我的NTT前面加了一个预处理的常数，所以前面FFT和FWT与NTT运行时间的比值会有较大的误差。

可以发现在实际运行中，所有算法与NTT的相对常数大概还要再乘1.51.51.5左右，其中任意模数NTT更是比理论复杂度慢了2倍（与FFT相比）。但是随着nnn的增大，时间复杂度相同的情况下，绝大多数比值在逐渐减小，这也符合渐进复杂度的规律。

多项式常数大，慎用！