多项式总结

多项式

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)=kln 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(nlog^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;
}

套路

对于(xin[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(nlog^2n))的。