多项式总结
多项式
FFT
void FFT(Complex *P,int opt){
for (int i=0;i<n;++i) if (i<rev[i]) swap(P[i],P[rev[i]]);
for (int i=1;i<n;i<<=1)
for (int p=i<<1,j=0;j<n;j+=p)
for (int k=0;k<i;++k){
Complex W=w[n/i*k];W.im*=opt;
Complex X=P[j+k],Y=W*P[j+k+i];
P[j+k]=X+Y;P[j+k+i]=X-Y;
}
if (opt==-1) for (int i=0;i<n;++i) P[i].rl/=1.0*n;
}
复数重载
struct Complex{
double rl,im;
Complex(){rl=im=0;}
Complex(double a,double b){rl=a,im=b;}
Complex operator + (Complex b)
{return Complex(rl+b.rl,im+b.im);}
Complex operator - (Complex b)
{return Complex(rl-b.rl,im-b.im);}
Complex operator * (Complex b)
{return Complex(rl*b.rl-im*b.im,rl*b.im+im*b.rl);}
}w[N];
单位根预处理
for (int i=0;i<n;++i) w[i]=Complex(cos(Pi*i/n),sin(Pi*i/n));
NTT
void NTT(int *P,int opt,int n){
int len,l=0;
for (len=1;len<n;len<<=1) ++l;--l;
for (int i=0;i<len;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<l);
for (int i=0;i<len;++i) if (i<rev[i]) swap(P[i],P[rev[i]]);
for (int i=1;i<len;i<<=1){
int W=fastpow(3,(mod-1)/(i<<1));
if (opt==-1) W=fastpow(W,mod-2);
og[0]=1;for (int j=1;j<i;++j) og[j]=1ll*og[j-1]*W%mod;
for (int p=i<<1,j=0;j<len;j+=p)
for (int k=0;k<i;++k){
int X=P[j+k],Y=1ll*P[j+k+i]*og[k]%mod;
P[j+k]=(X+Y)%mod;P[j+k+i]=(X-Y+mod)%mod;
}
}
if (opt==-1)
for (int i=0,Inv=fastpow(len,mod-2);i<len;++i)
P[i]=1ll*P[i]*Inv%mod;
}
MTT
求\(F(x)\)与\(G(x)\)在任意模数下的卷积。
为什么不能直接\(FFT\)乘然后再取模？因为直接乘结果会爆long long。
考虑拆系数。设一个常数\(M\)，把\(F(x)\)和\(G(x)\)拆成：\(A(x)=\frac{F(x)}M\)，\(B(x)=F(x)\%M\)，\(C(x)=\frac{G(x)}M\)，\(D(x)=G(x)\%M\)。
这样答案就变成了：
直接\(FFT\)就不会爆long long了。\(M\)取\(\sqrt {mod}\)最优。
这种方法一共要做7次\(DFT\)。
for (int i=0;i<=n;++i){
int x=gi()%mod;
a[i]=Complex(x/qmod,0);b[i]=Complex(x%qmod,0);
}
for (int i=0;i<=m;++i){
int x=gi()%mod;
c[i]=Complex(x/qmod,0);d[i]=Complex(x%qmod,0);
}
FFT(a,1);FFT(b,1);FFT(c,1);FFT(d,1);
for (int i=0;i<n;++i){
s1[i]=a[i]*c[i];
s2[i]=a[i]*d[i]+b[i]*c[i];
s3[i]=b[i]*d[i];
}
FFT(s1,-1);FFT(s2,-1);FFT(s3,-1);
for (int i=0;i<=m;++i){
int ans=(ll)(s1[i].rl+0.5)%mod*qmod%mod*qmod%mod;
(ans+=(ll)(s2[i].rl+0.5)%mod*qmod%mod)%=mod;
(ans+=(ll)(s3[i].rl+0.5)%mod)%=mod;
printf("%d ",ans);
}
多项式运算
假设下面都是给出\(A(x)\)要你求\(B(x)\)
所以可以把\(A(x)\)看作是一个常数而\(B(x)\)是函数的自变量
牛顿迭代
很多东西就只要套式子就可以了。注意这里的求导是对\(B(x)\)求导而不是对\(x\)求导。
多项式求导
有\((x^a)'=ax^{a-1}\)，而且导数是满足线性性的。
所以
void Dao(int *a,int *b,int len){
for (int i=1;i<len;++i) b[i-1]=1ll*i*a[i]%mod;
b[len]=b[len-1]=0;
}
多项式求积分
\(\int x^a=\frac{1}{a+1}x^{a+1}\)，积分同样满足线性性。
void Jifen(int *a,int *b,int len){
for (int i=1;i<len;++i) b[i]=1ll*a[i-1]*inv[i]%mod;
b[0]=0;
}
多项式求逆：
int A[_],B[_];
void GetInv(int *a,int *b,int len){
if (len==1) {b[0]=fastpow(a[0],mod-2);return;}
GetInv(a,b,len>>1);
for (int i=0;i<len;++i) A[i]=a[i],B[i]=b[i];
NTT(A,1,len<<1);NTT(B,1,len<<1);
for (int i=0;i<(len<<1);++i) A[i]=1ll*A[i]*B[i]%mod*B[i]%mod;
NTT(A,-1,len<<1);
for (int i=0;i<len;++i) b[i]=((b[i]+b[i])%mod-A[i]+mod)%mod;
for (int i=0;i<(len<<1);++i) A[i]=B[i]=0;
}
多项式开方
int C[_],D[_],inv2=fastpow(2,mod-2);
void GetSqrt(int *a,int *b,int len){
if (len==1) {b[0]=sqrt(a[0]);return;}
GetSqrt(a,b,len>>1);
for (int i=0;i<len;++i) C[i]=a[i];
GetInv(b,D,len);
NTT(C,1,len<<1);NTT(D,1,len<<1);
for (int i=0;i<(len<<1);++i) D[i]=1ll*D[i]*C[i]%mod;
NTT(D,-1,len<<1);
for (int i=0;i<len;++i) b[i]=1ll*(b[i]+D[i])%mod*inv2%mod;
for (int i=0;i<(len<<1);++i) C[i]=D[i]=0;
}
多项式求\(\ln\)：
相当于先对\(A(x)\)求个导求个逆然后乘起来再积分
void Getln(int *a,int *b,int len){
int A[_],B[_];
memset(A,0,sizeof(A));memset(B,0,sizeof(B));
Dao(a,A,len);GetInv(a,B,len);
NTT(A,1,len<<1);NTT(B,1,len<<1);
for (int i=0;i<(len<<1);++i) A[i]=1ll*A[i]*B[i]%mod;
NTT(A,-1,len<<1);
Jifen(A,b,len);
}
注意多项式求\(\ln\)要保证\(A(x)\)常数项为\(1\)，求完后默认常数项为\(0\)。
多项式\(\exp\)：
int D[_],E[_];
void GetExp(int *a,int *b,int len){
if (len==1) {b[0]=1;return;}
GetExp(a,b,len>>1);
for (int i=0;i<len;++i) D[i]=b[i];
Getln(b,E,len);
for (int i=0;i<len;++i) E[i]=(mod-E[i]+a[i])%mod;E[0]=(E[0]+1)%mod;
NTT(D,1,len<<1);NTT(E,1,len<<1);
for (int i=0;i<(len<<1);++i) D[i]=1ll*D[i]*E[i]%mod;
NTT(D,-1,len<<1);
for (int i=0;i<len;++i) b[i]=D[i];
for (int i=0;i<(len<<1);++i) D[i]=E[i]=0;
}
注意多项式求\(\exp\)要保证\(A(x)\)常数项为\(0\)，求完后默认常数项为\(1\)。
多项式快速幂
相当于先取对数，然后每个系数乘以下，再做个\(\exp\)。
void GetPow(int *a,int *b,int len){
int F[_];memset(F,0,sizeof(F));
Getln(a,F,len);
for (int i=0;i<len;++i) F[i]=1ll*F[i]*k%mod;
GetExp(F,b,len);
}
分治FFT
不是\(CDQ\)的那套理论，就是单纯计算\(\prod_{i=1}^{n}(1+a_ix)\)
理论上讲起来挺容易的，复杂度\(O(n\log^2n)\)。
int tmp[50][_],Stack[50],top;
void solve(int *P,int *q,int l,int r){
if (l==r){P[0]=1;P[1]=q[l];return;}
int mid=l+r>>1,ls=Stack[top--];
solve(tmp[ls],q,l,mid);
int rs=Stack[top--];
solve(tmp[rs],q,mid+1,r);
int len=1;
while (len<=r-l+1) len<<=1;
NTT(tmp[ls],1,len);NTT(tmp[rs],1,len);
for (int i=0;i<len;++i) P[i]=1ll*tmp[ls][i]*tmp[rs][i]%mod;
NTT(P,-1,len);
Stack[++top]=ls;Stack[++top]=rs;
for (int i=0;i<len;++i) tmp[ls][i]=tmp[rs][i]=0;
}
套路
对于\(x\in[1,m]\)计算\(\sum_{i=1}^{n}a_i^x\)。
先用分治\(FFT\)计算一下上面的那个东西也就是\(f(x)=\prod_{i=1}^{n}(1+a_ix)\)。
取对数\(\ln f(x)=\sum_{i=1}^{n}\ln(1+a_ix)\)
对这个东西求个导。
\((\ln f(x))'=\sum_{i=1}^{n}(\ln(1+a_ix))'=\sum_{i=1}^{n}\frac{a_i}{1+a_ix}\)。
这是后面是一个无限项等比数列求和的形式。
\((\ln f(x))'=\sum_{i=1}^{n}\sum_{j=0}^{\inf}(-1)^ja_i^{j+1}x^j=\sum_{j=0}^{\inf}(-1)^j(\sum_{i=1}^{n}a_i^{j+1})x_j\)。
所以求出\(f(x)\)以后求\(\ln\)求导再把有的项取反就行了。注意这里求的是\(\sum_{i=1}^{n}a_i^{j+1}\)，如果你要求\(\sum_{i=1}^n a_i^0\)的话，那不就是\(n\)吗。
这种做法是\(O(n\log^2n)\)的。
