多项式总结

多项式

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\)。
这样答案就变成了：

\[M^2A(x)C(x)+M(A(x)D(x)+B(x)C(x))+B(x)D(x) \]

直接\(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_{t+1}(x)=B_t(x)-\frac{F(B_t(x))}{F'(B_t(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;
}

多项式求逆：

\[A(x)B(x)=1\\A(x)B(x)-1=0\\B_{t+1}(x)=B_t(x)-\frac{A(x)B_t(x)-1}{A(x)}\\B_{t+1}(x)=B_t(x)-B_t(x)(A(x)B_t(x)-1)\\B_{t+1}(x)=2B_t(x)-A(x)B_t^2(x) \]

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

多项式开方

\[A(x)=B(x)B(x)\\B_{t+1}(x)=B_t(x)-\frac{B_t(x)B_t(x)-A(x)}{2B_t(x)}\\B_{t+1}(x)=\frac 12(\frac{A(x)}{B(x)}+B(x)) \]

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\)：

\[\ln A(x)=B(x)\\\frac{A'(x)}{A(x)}=B'(x)\\ \]

相当于先对\(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\)：

\[e^{A(x)}=B(x)\\A(x)=\ln B(x)\\\ln B(x)-A(x)=0\\B_{t+1}(x)=B_t(x)-\frac{\ln B(x)-A(x)}{\frac{1}{B(x)}}\\B_{t+1}(x)=B_t(x)(1-\ln B(x)+A(x)) \]

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\)。

多项式快速幂

\[B(x)=A^k(x)\\\ln B(x)=k\ln A(x) \]

相当于先取对数，然后每个系数乘以下，再做个\(\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)\)的。