多项式求逆 学习总结

感觉蒟蒻现在学多项式求逆貌似有些晚了

不过真的很有意思了(然而省选的时候自己还在玩泥巴什么也不会

多项式求逆是基于倍增的

假设我们知道

h(x)f(x)=1(mod x^n)

移项得

(h(x)*f(x)-1)=0(mod x^n)

两边同时求平方得

h(x)^2*f(x)^2 - 2*h(x)*f(x) +1 = 0 (mod x^2n)

设g(x)*f(x)=1(mod x^2n)

两边同时乘以g(x)可以得

h(x)^2*f(x) -2*h(x) + g(x) =0 （mod x^2n)

我们移项可以得到

g(x) = h(x) *(2- f(x)*h(x)) 

倍增即可，时间复杂度O(nlogn)

首先是picks的blog里提到的伯努利数

由于并没有找到相应的题目，于是就在cojs上自己出了一发

疯狂的求和问题 80分

我们知道伯努利数的生成函数是x/(e^x-1)

我们把下面做泰勒展开之后进行多项式求逆即可

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#define G 3
using namespace std;

typedef long long LL;
const int mod=998244353;
const int maxn=500010;
int n,k,N,len,ans;
int jc[maxn],inv[maxn];
int rev[maxn];
int A[maxn],B[maxn],C[maxn];

void read(int &num){
	num=0;char ch=getchar();
	while(ch<'!')ch=getchar();
	while(ch>='0'&&ch<='9'){
		num=(10LL*num+ch-'0')%mod;
		ch=getchar();
	}return;
}
int pow_mod(int v,int p){
	int tmp=1;
	while(p){
		if(p&1)tmp=1LL*tmp*v%mod;
		v=1LL*v*v%mod;p>>=1;
	}return tmp;
}
void FFT(int *A,int n,int flag){
	for(len=0;(1<<len)<n;++len);
	for(int i=0;i<n;++i){
		rev[i]=rev[i>>1]>>1|((i&1)<<(len-1));
		C[i]=A[rev[i]];
	}
	for(int i=0;i<n;++i)A[i]=C[i];
	for(int k=2;k<=n;k<<=1){
		int mk=(k>>1);
		int wn=pow_mod(G,flag==1?(mod-1)/k:(mod-1)-(mod-1)/k);
		for(int i=0;i<n;i+=k){
			int w=1;
			for(int j=0;j<mk;++j){
				int a=i+j,b=i+j+mk;
				int x=A[a],y=1LL*A[b]*w%mod;
				A[a]=x+y;if(A[a]>=mod)A[a]-=mod;
				A[b]=x-y;if(A[b]<0)A[b]+=mod;
				w=1LL*w*wn%mod;
			}
		}
	}
	if(flag==-1){
		int c=pow_mod(n,mod-2);
		for(int i=0;i<n;++i)A[i]=1LL*A[i]*c%mod;
	}return;
}
void Get_inv(int n){
	if(n==1){
		B[0]=pow_mod(A[0],mod-2);
		return;
	}
	int now=(n>>1);Get_inv(now);
	static int tmp[maxn];
	for(int i=0;i<n;++i)tmp[i]=A[i];
	now=(n<<1);
	FFT(B,now,1);FFT(tmp,now,1);
	for(int i=0;i<now;++i)B[i]=1LL*B[i]*(2LL-1LL*tmp[i]*B[i]%mod+mod)%mod;
	FFT(B,now,-1);
	memset(B+n,0,sizeof(int)*n);
}
int Get_C(int n,int m){return 1LL*jc[n]*inv[m]%mod*inv[n-m]%mod;}

int main(){
	freopen("Crazy_Sum.in","r",stdin);
	freopen("Crazy_Sum.out","w",stdout);
	read(n);scanf("%d",&k);
	for(N=1;N<=k;N<<=1);
	jc[0]=1;
	for(int i=1;i<=N;++i)jc[i]=1LL*jc[i-1]*i%mod;
	inv[N]=pow_mod(jc[N],mod-2);
	for(int i=N-1;i>=0;--i)inv[i]=1LL*inv[i+1]*(i+1)%mod;
	for(int i=0;i<N;++i)A[i]=inv[i+1];
	Get_inv(N);
	for(int i=0;i<N;++i)B[i]=1LL*B[i]*jc[i]%mod;
	int now=n+1;
	for(int i=1;i<=k+1;++i){
		ans=ans+1LL*Get_C(k+1,i)*B[k+1-i]%mod*now%mod;
		if(ans>=mod)ans-=mod;
		now=1LL*now*(n+1)%mod;
	}ans=1LL*ans*pow_mod(k+1,mod-2)%mod;
	printf("%d
",ans);
	return 0;
}

BZOJ 3456

很经典的题目啦，设fi表示i个点的联通图的个数，考虑1号点所在联通块的大小我们有

fn=2^C(n,2)-sigma(C(n-1,i-1)*fi*2^C(n-i,2))

很容易发现这是个CDQ+FFT的式子，直接上就可以了

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#define G 3
using namespace std;
 
const int mod=1004535809;
const int maxn=300010;
int n,N,len;
int jc[maxn],inv[maxn],rev[maxn];
int h[maxn],f[maxn];
int A[maxn],B[maxn],C[maxn];
 
int pow_mod(int v,int p){
    int tmp=1;
    while(p){
        if(p&1)tmp=1LL*tmp*v%mod;
        v=1LL*v*v%mod;p>>=1;
    }return tmp;
}
void FFT(int *A,int n,int flag){
    for(int i=0;i<n;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
    for(int k=2;k<=n;k<<=1){
        int mk=(k>>1);
        int wn=pow_mod(G,flag==1?(mod-1)/k:(mod-1)-(mod-1)/k);
        for(int i=0;i<n;i+=k){
            int w=1;
            for(int j=0;j<mk;++j){
                int a=i+j,b=i+j+mk;
                int x=A[a],y=1LL*A[b]*w%mod;
                A[a]=x+y;if(A[a]>=mod)A[a]-=mod;
                A[b]=x-y;if(A[b]<0)A[b]+=mod;
                w=1LL*w*wn%mod;
            }
        }
    }
    if(flag==-1){
        int c=pow_mod(n,mod-2);
        for(int i=0;i<n;++i)A[i]=1LL*A[i]*c%mod;
    }return;
}
void Solve(int L,int R){
    if(L==R)return;
    int mid=(L+R)>>1;
    Solve(L,mid);
    for(N=1,len=0;N<=(R-L+1);N<<=1,len++);
    for(int i=0;i<N;++i)A[i]=0,B[i]=h[i];
    for(int i=L;i<=mid;++i)A[i-L]=1LL*f[i]*inv[i-1]%mod;
    if(R-L<=500){
        int lim=R-L+1;
        for(int i=0;i<lim;++i){
            C[i]=0;
            for(int j=0;j<=i;++j){
                C[i]=C[i]+1LL*A[j]*B[i-j]%mod;
                if(C[i]>=mod)C[i]-=mod;
            }
        }
    }else{
        for(int i=0;i<N;++i)rev[i]=rev[i>>1]>>1|((i&1)<<(len-1));
        FFT(A,N,1);FFT(B,N,1);
        for(int i=0;i<N;++i)C[i]=1LL*A[i]*B[i]%mod;
        FFT(C,N,-1);
    }
    for(int i=mid+1;i<=R;++i){
        f[i]=f[i]-1LL*jc[i-1]*C[i-L]%mod;
        if(f[i]<0)f[i]+=mod;
    }
    Solve(mid+1,R);
}
 
int main(){
    scanf("%d",&n);
    jc[0]=1;
    for(int i=1;i<=n;++i)jc[i]=1LL*jc[i-1]*i%mod;
    inv[n]=pow_mod(jc[n],mod-2);
    for(int i=n-1;i>=0;--i)inv[i]=1LL*inv[i+1]*(i+1)%mod;
    h[1]=1;f[1]=1;
    for(int i=2;i<=n;++i){
        int tmp=(1LL*i*(i-1)/2)%(mod-1);
        h[i]=pow_mod(2,tmp);f[i]=h[i];
        h[i]=1LL*h[i]*inv[i]%mod;
    }
    Solve(1,n);
    printf("%d
",f[n]);
    return 0;
}

但是这样我们是两个log的

我们注意到可以把上面的式子等价转化成

2^C(n,2)/(n-1)!=sigma( 2^C(n-i,2)/(n-i)! *fi/(i-1)! )

可以构造出h=g*f的形式

之后我们已知h和g，求f，求h*g^(-1)即可,也就是多项式求逆

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<algorithm>
#define G 3
using namespace std;
 
typedef long long LL;
const int maxn=500010;
const int mod=1004535809;
int n,N,len,t;
int jc[maxn],inv[maxn],rev[maxn];
int A[maxn],B[maxn];
 
int pow_mod(int v,int p){
    int tmp=1;
    while(p){
        if(p&1)tmp=1LL*tmp*v%mod;
        v=1LL*v*v%mod;p>>=1;
    }return tmp;
}
void FFT(int *A,int n,int flag){
    for(int i=0;i<n;++i){
        if(i<rev[i]){t=A[i];A[i]=A[rev[i]];A[rev[i]]=t;}
    }
    for(int k=2;k<=n;k<<=1){
        int mk=(k>>1);
        int wn=pow_mod(G,flag==1?(mod-1)/k:(mod-1)-(mod-1)/k);
        for(int i=0;i<n;i+=k){
            int w=1;
            for(int j=0;j<mk;++j){
                int a=i+j,b=i+j+mk;
                int x=A[a],y=1LL*A[b]*w%mod;
                A[a]=x+y;if(A[a]>=mod)A[a]-=mod;
                A[b]=x-y;if(A[b]<0)A[b]+=mod;
                w=1LL*w*wn%mod;
            }
        }
    }
    if(flag==-1){
        int c=pow_mod(n,mod-2);
        for(int i=0;i<n;++i)A[i]=1LL*A[i]*c%mod;
    }return;
}
void Get_inv(int n){
    if(n==1){
        B[0]=pow_mod(A[0],mod-2);
        return;
    }
    Get_inv(n>>1);
    int now=(n<<1);
    for(len=0;(1<<len)<now;++len);
    for(int i=0;i<now;++i)rev[i]=rev[i>>1]>>1|((i&1)<<(len-1));
    static int tmp[maxn];
    for(int i=0;i<n;++i)tmp[i]=A[i];
    FFT(tmp,now,1);FFT(B,now,1);
    for(int i=0;i<now;++i)B[i]=1LL*B[i]*(2LL-1LL*tmp[i]*B[i]%mod+mod)%mod;
    FFT(B,now,-1);
    memset(B+n,0,sizeof(int)*n);
}
 
int main(){
    scanf("%d",&n);
    for(N=1;N<=n;N<<=1);
    jc[0]=1;
    for(int i=1;i<N;++i)jc[i]=1LL*jc[i-1]*i%mod;
    inv[N-1]=pow_mod(jc[N-1],mod-2);
    for(int i=N-2;i>=0;--i)inv[i]=1LL*inv[i+1]*(i+1)%mod;
    for(int i=0;i<N;++i){
        int tmp=(1LL*i*(i-1)/2)%(mod-1);
        A[i]=1LL*pow_mod(2,tmp)*inv[i]%mod;
    }
    Get_inv(N);
    memset(B+n,0,sizeof(int)*n);
    for(len=0;(1<<len)<N;++len);
    for(int i=0;i<N;++i)rev[i]=rev[i>>1]>>1|((i&1)<<(len-1));
    for(int i=0;i<N;++i){
        int tmp=(1LL*i*(i-1)/2)%(mod-1);
        A[i]=1LL*pow_mod(2,tmp)*inv[i-1]%mod;
    }
    for(len=0;(1<<len)<N;++len);
    for(int i=0;i<N;++i)rev[i]=rev[i>>1]>>1|((i&1)<<(len-1));
    FFT(A,N,1);FFT(B,N,1);
    for(int i=0;i<N;++i)A[i]=1LL*A[i]*B[i]%mod;
    FFT(A,N,-1);
    printf("%d
",1LL*A[n]*jc[n-1]%mod);
    return 0;
}

cojs 疯狂的欧拉图

搞完了城市规划之后自己把原来的考试题扩大了数据范围出到了cojs上(原来n^2可过)

式子还是一样的推导，之后也可以转换成多项式求逆

但是由于模数不兹瓷，所以我们要写模任意质数的FFT QAQ

CDQ+FFT

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#define G 3
using namespace std;
 
const int maxn=400010;
typedef long long LL;
int n,N,len;
const int mod=1e9+7;
const int M=30000;
const long double pi=acos(-1.0);
int jc[maxn],inv[maxn];
int h[maxn],f[maxn];
int rev[maxn];
int a[maxn],b[maxn],c[maxn];
int a0[maxn],b0[maxn],a1[maxn],b1[maxn];
struct cpx{
    long double r,i;
    cpx(long double r=0,long double i=0):r(r),i(i){}
}A[maxn],B[maxn],C[maxn];
cpx operator +(const cpx &A,const cpx &B){return cpx(A.r+B.r,A.i+B.i);}
cpx operator -(const cpx &A,const cpx &B){return cpx(A.r-B.r,A.i-B.i);}
cpx operator *(const cpx &A,const cpx &B){return cpx(A.r*B.r-A.i*B.i,A.r*B.i+A.i*B.r);}
void FFT(cpx *A,int n,int type){
    for(int i=0;i<n;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
    for(int k=0;(1<<k)<n;++k){
        int m=(1<<k),m2=(m<<1);
        long double o=2*pi/m2*type;
        cpx wn(cos(o),sin(o));
        for(int i=0;i<n;i+=m2){
            cpx w(1,0);
            for(int j=0;j<m;++j){
                cpx x=A[i+j],y=A[i+j+m]*w;
                A[i+j]=x+y;A[i+j+m]=x-y;
                w=w*wn;
            }
        }
    }
    if(type==-1)for(int i=0;i<n;++i)A[i].r/=n;
}
void mul(int *a,int *b,int *c){
    for(int i=0;i<N;++i)A[i]=cpx(a[i],0),B[i]=cpx(b[i],0);
    FFT(A,N,1);FFT(B,N,1);
    for(int i=0;i<N;++i)A[i]=A[i]*B[i];
    FFT(A,N,-1);
    for(int i=0;i<N;++i)c[i]=((LL)(A[i].r+0.5))%mod;
}
void mul_mod(int *a,int *b,int *c){
    for(int i=0;i<N;++i)a0[i]=a[i]/M,b0[i]=b[i]/M;
    mul(a0,b0,a0);
    for(int i=0;i<N;++i){
        c[i]=1LL*a0[i]*M*M%mod;
        a1[i]=a[i]%M;b1[i]=b[i]%M;
    }
    mul(a1,b1,a1);
    for(int i=0;i<N;++i){
        c[i]=c[i]+a1[i];
        if(c[i]>=mod)c[i]-=mod;
        a1[i]=a1[i]+a0[i];
        if(a1[i]>=mod)a1[i]-=mod;
        a0[i]=a[i]/M+a[i]%M;
        b0[i]=b[i]/M+b[i]%M;
    }
    mul(a0,b0,a0);
    for(int i=0;i<N;++i){
        c[i]=c[i]+1LL*M*(a0[i]-a1[i]+mod)%mod;
        if(c[i]>=mod)c[i]-=mod;
    }return;
}
LL pow_mod(LL v,int p){
    LL tmp=1;
    while(p){
        if(p&1)tmp=tmp*v%mod;
        v=v*v%mod;p>>=1;
    }return tmp;
}
void Solve(int L,int R){
    if(L==R)return;
    int mid=(L+R)>>1;
    Solve(L,mid);
    for(N=1,len=0;N<(R-L+1);N<<=1,len++);
    for(int i=0;i<N;++i)a[i]=b[i]=0;
    for(int i=L;i<=mid;++i)a[i-L]=1LL*f[i]*inv[i-1]%mod;
    for(int i=0;i<N;++i)b[i]=h[i];
    if(R-L<=1000){
    	int lim=R-L+1;
    	for(int i=0;i<lim;++i){
    		c[i]=0;
    		for(int j=0;j<=i;++j){
    			c[i]=c[i]+1LL*a[j]*b[i-j]%mod;
    			if(c[i]>=mod)c[i]-=mod;
    		}
    	}
    }else{
    	for(int i=0;i<N;++i)rev[i]=rev[i>>1]>>1|((i&1)<<(len-1));
    	mul_mod(a,b,c);
    }
    for(int i=mid+1;i<=R;++i){
        f[i]=f[i]-1LL*jc[i-1]*c[i-L]%mod;
        if(f[i]<0)f[i]+=mod;
    }
    Solve(mid+1,R);
}
 
int main(){
	freopen("Crazy_Graph.in","r",stdin);
	freopen("Crazy_Graph.out","w",stdout);
    scanf("%d",&n);
    jc[0]=1;
    for(int i=1;i<=n;++i)jc[i]=1LL*jc[i-1]*i%mod;
    inv[n]=pow_mod(jc[n],mod-2);
    for(int i=n-1;i>=0;--i)inv[i]=1LL*inv[i+1]*(i+1)%mod;
    h[1]=1;f[1]=1;
    for(int i=2;i<=n;++i){
        h[i]=pow_mod(2LL,(1LL*(i-1)*(i-2)/2)%(mod-1));
        f[i]=h[i];
        h[i]=1LL*h[i]*inv[i]%mod;
    }
    Solve(1,n);
    f[n]=f[n]+1LL*n*(n-1)*pow_mod(2LL,mod-2)%mod*f[n]%mod;
    if(f[n]>=mod)f[n]-=mod;
    printf("%d
",f[n]);
    return 0;
}

多项式求逆

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<cmath>
using namespace std;

typedef long long LL;
const int mod=1e9+7;
const int M=32000;
const int maxn=400010;
const long double pi=acos(-1.0);
int n,len,N;
int jc[maxn],inv[maxn],rev[maxn];
int a[maxn],b[maxn],c[maxn];
int a0[maxn],b0[maxn],a1[maxn],b1[maxn];
struct cpx{
	long double r,i;
	cpx(long double r=0,long double i=0):r(r),i(i){}
}A[maxn],B[maxn];
cpx operator +(const cpx &A,const cpx &B){return cpx(A.r+B.r,A.i+B.i);}
cpx operator -(const cpx &A,const cpx &B){return cpx(A.r-B.r,A.i-B.i);}
cpx operator *(const cpx &A,const cpx &B){return cpx(A.r*B.r-A.i*B.i,A.r*B.i+A.i*B.r);}

int pow_mod(int v,int p){
	int tmp=1;
	while(p){
		if(p&1)tmp=1LL*tmp*v%mod;
		v=1LL*v*v%mod;p>>=1;
	}return tmp;
}
void FFT(cpx *A,int n,int flag){
	for(int i=0;i<n;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
	for(int k=2;k<=n;k<<=1){
		int mk=(k>>1);
		long double o=2*pi/k*flag;
		cpx wn(cos(o),sin(o));
		for(int i=0;i<n;i+=k){
			cpx w(1,0);
			for(int j=0;j<mk;++j){
				int a=i+j,b=i+j+mk;
				cpx x=A[a],y=A[b]*w;
				A[a]=x+y;A[b]=x-y;
				w=w*wn;
			}
		}
	}
	if(flag==-1){for(int i=0;i<n;++i)A[i].r/=n;}
}
void mul(int *a,int *b,int *c,int N){
	for(int i=0;i<N;++i)A[i]=cpx(a[i],0),B[i]=cpx(b[i],0);
	FFT(A,N,1);FFT(B,N,1);
	for(int i=0;i<N;++i)A[i]=A[i]*B[i];
	FFT(A,N,-1);
	for(int i=0;i<N;++i)c[i]=(LL)(A[i].r+0.5)%mod;
}
void mul_mod(int *A,int *B,int *C,int N){
	for(int i=0;i<N;++i)a0[i]=A[i]%M,b0[i]=B[i]%M;
	mul(a0,b0,a0,N);
	for(int i=0;i<N;++i)C[i]=a0[i],a1[i]=A[i]/M,b1[i]=B[i]/M;
	mul(a1,b1,a1,N);
	for(int i=0;i<N;++i){
		C[i]=C[i]+1LL*M*M*a1[i]%mod;
		if(C[i]>=mod)C[i]-=mod;
		a1[i]=a1[i]+a0[i];
		if(a1[i]>=mod)a1[i]-=mod;
		a0[i]=A[i]/M+A[i]%M;
		b0[i]=B[i]/M+B[i]%M;
	}
	mul(a0,b0,a0,N);
	for(int i=0;i<N;++i){
		C[i]=C[i]+1LL*M*(a0[i]-a1[i]+mod)%mod;
		if(C[i]>=mod)C[i]-=mod;
	}return;
}
void Get_inv(int n){
	if(n==1){
		b[0]=pow_mod(a[0],mod-2);
		return;
	}
	Get_inv(n>>1);
	int now=(n<<1);
	for(len=0;(1<<len)<now;++len);
	for(int i=0;i<now;++i)rev[i]=rev[i>>1]>>1|((i&1)<<(len-1));
	static int tmp[maxn];
	for(int i=0;i<n;++i)tmp[i]=a[i];
	mul_mod(b,tmp,c,now);
	c[0]=2-c[0];
	if(c[0]<0)c[0]+=mod;
	for(int i=1;i<now;++i)c[i]=mod-c[i];
	mul_mod(b,c,tmp,now);
	for(int i=0;i<n;++i)b[i]=tmp[i];
	memset(b+n,0,sizeof(int)*n);
}

int main(){
	freopen("Crazy_Graph.in","r",stdin);
	freopen("Crazy_Graph.out","w",stdout);
	scanf("%d",&n);
	for(N=1;N<=n;N<<=1);
	jc[0]=1;
	for(int i=1;i<N;++i)jc[i]=1LL*jc[i-1]*i%mod;
	inv[N-1]=pow_mod(jc[N-1],mod-2);
	for(int i=N-2;i>=0;--i)inv[i]=1LL*inv[i+1]*(i+1)%mod;
	a[0]=1;
	for(int i=1;i<N;++i){
		int tmp=(1LL*(i-1)*(i-2)/2)%(mod-1);
		a[i]=1LL*pow_mod(2,tmp)*inv[i]%mod;
	}
	Get_inv(N);
	memset(b+n,0,sizeof(int)*n);
	for(len=0;(1<<len)<N;++len);
	for(int i=0;i<N;++i)rev[i]=rev[i>>1]>>1|((i&1)<<(len-1));
	a[0]=1;
	for(int i=1;i<N;++i){
		int tmp=(1LL*(i-1)*(i-2)/2)%(mod-1);
		a[i]=1LL*pow_mod(2,tmp)*inv[i-1]%mod;
	}
	mul_mod(a,b,c,N);
	int ans=1LL*jc[n-1]*c[n]%mod;
	ans=ans+1LL*n*(n-1)*pow_mod(2LL,mod-2)%mod*ans%mod;
	if(ans>=mod)ans-=mod;
	printf("%d
",ans);
	return 0;
}

cojs 异化多肽

貌似是多项式求逆的裸题，Nescafe的题目

构造生成函数f，不难发现我们要求的是f^1+f^2+……

之后求和得 1/(1-f)

多项式求逆即可

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cstdlib>
#include<algorithm>
#define G 5
using namespace std;

typedef long long LL;
const int maxn=300010;
const int mod=1005060097;
int n,m,k,N,L;
int A[maxn],B[maxn];
int rev[maxn];

void read(int &num){
	num=0;char ch=getchar();
	while(ch<'!')ch=getchar();
	while(ch>='0'&&ch<='9')num=num*10+ch-'0',ch=getchar();
}
int pow_mod(int v,int p){
	int tmp=1;
	while(p){
		if(p&1)tmp=1LL*tmp*v%mod;
		v=1LL*v*v%mod;p>>=1;
	}return tmp;
}
void FFT(int *A,int n,int flag){
	for(L=0;(1<<L)<n;++L);
	for(int i=0;i<n;++i)rev[i]=rev[i>>1]>>1|((i&1)<<(L-1));
	for(int i=0;i<n;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
	for(int k=2;k<=n;k<<=1){
		int mk=(k>>1);
		int wn=pow_mod(G,flag==1?(mod-1)/k:(mod-1)-(mod-1)/k);
		for(int i=0;i<n;i+=k){
			int w=1;
			for(int j=0;j<mk;++j){
				int x=A[i+j],y=1LL*A[i+j+mk]*w%mod;
				A[i+j]=x+y;if(A[i+j]>=mod)A[i+j]-=mod;
				A[i+j+mk]=x-y;if(A[i+j+mk]<0)A[i+j+mk]+=mod;
				w=1LL*w*wn%mod;
			}
		}
	}
	if(flag==-1){
		int inv=pow_mod(n,mod-2);
		for(int i=0;i<n;++i)A[i]=1LL*A[i]*inv%mod;
	}return;
}
void Get_inv(int n){
	if(n==1){
		B[0]=pow_mod(A[0],mod-2);
		return;
	}
	int now=(n>>1);
	Get_inv(now);
	static int tmp[maxn];
	for(int i=0;i<now;++i)tmp[i]=A[i];
	FFT(tmp,n,1);FFT(B,n,1);
	for(int i=0;i<n;++i)B[i]=1LL*B[i]*(2-1LL*tmp[i]*B[i]%mod+mod)%mod;
	FFT(B,n,-1);
	memset(B+now,0,sizeof(int)*now);
}

int main(){
	freopen("polypeptide.in","r",stdin);
	freopen("polypeptide.out","w",stdout);
	read(n);read(m);
	for(int i=1;i<=m;++i){
		read(k);
		if(k<=n)A[k]++;
	}A[0]--;if(A[0]<0)A[0]+=mod;
	for(N=1;N<=n;N<<=1);N<<=1;
	Get_inv(N);
	int ans=-B[n];
	ans=(ans%mod+mod)%mod;
	printf("%d
",ans);
	return 0;
}

HEOI 2016 求和

QAQ 当前考场上没推出式子来

其实当时并不是很会FFT，现在也不是很会

我们不难发现如果没有2^j的话求的是n个数划分成若干个有序集合的方案数

设fn为这个方案数，考虑枚举第一个集合的大小

我们有fn=sigma( C(n,i) * f(n-i) )

考虑2^j实际上是每个集合都要*2

那么我们的递推式只要改成 fn=sigma( C(n,i) * f(n-i) *2)就可以了

然后上CDQ+FFT

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<algorithm>
#define G 3
using namespace std;

typedef long long LL;
const int maxn=500010;
const int mod=998244353;
int n,N,len,ans;
int jc[maxn],inv[maxn],rev[maxn];
int f[maxn],A[maxn],B[maxn],C[maxn];

int pow_mod(int v,int p){
	int tmp=1;
	while(p){
		if(p&1)tmp=1LL*tmp*v%mod;
		v=1LL*v*v%mod;p>>=1;
	}return tmp;
}
void FFT(int *A,int n,int flag){
	for(int i=0;i<n;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
	for(int k=2;k<=n;k<<=1){
		int mk=(k>>1);
		int wn=pow_mod(G,flag==1?(mod-1)/k:(mod-1)-(mod-1)/k);
		for(int i=0;i<n;i+=k){
			int w=1;
			for(int j=0;j<mk;++j){
				int a=i+j,b=i+j+mk;
				int x=A[a],y=1LL*A[b]*w%mod;
				A[a]=x+y;if(A[a]>=mod)A[a]-=mod;
				A[b]=x-y;if(A[b]<0)A[b]+=mod;
				w=1LL*w*wn%mod;
			}
		}
	}
	if(flag==-1){
		int c=pow_mod(n,mod-2);
		for(int i=0;i<n;++i)A[i]=1LL*A[i]*c%mod;
	}return;
}
void Solve(int L,int R){
	if(L==R)return;
	int mid=(L+R)>>1;
	Solve(L,mid);
	for(N=1,len=0;N<=(R-L+1);N<<=1,len++);
	for(int i=0;i<N;++i)A[i]=0,B[i]=inv[i];
	for(int i=L;i<=mid;++i)A[i-L]=f[i];
	if(R-L<=500){
		int lim=R-L+1;
		for(int i=0;i<lim;++i){
			C[i]=0;
			for(int j=0;j<=i;++j){
				C[i]=C[i]+1LL*A[j]*B[i-j]%mod;
				if(C[i]>=mod)C[i]-=mod;
			}
		}
	}else{
		for(int i=0;i<N;++i)rev[i]=rev[i>>1]>>1|((i&1)<<(len-1));
		FFT(A,N,1);FFT(B,N,1);
		for(int i=0;i<N;++i)C[i]=1LL*A[i]*B[i]%mod;
		FFT(C,N,-1);
	}
	for(int i=mid+1;i<=R;++i){
		f[i]=f[i]+2LL*C[i-L]%mod;
		if(f[i]>=mod)f[i]-=mod;
	}
	Solve(mid+1,R);
}

int main(){
	freopen("heoi2016_sum.in","r",stdin);
	freopen("heoi2016_sum.out","w",stdout);
	scanf("%d",&n);
	jc[0]=1;
	for(int i=1;i<=n;++i)jc[i]=1LL*jc[i-1]*i%mod;
	inv[n]=pow_mod(jc[n],mod-2);f[0]=1;
	for(int i=n-1;i>=0;--i)inv[i]=1LL*inv[i+1]*(i+1)%mod;
	Solve(0,n);
	for(int i=0;i<=n;++i){
		ans=ans+1LL*jc[i]*f[i]%mod;
		if(ans>=mod)ans-=mod;
	}printf("%d
",ans);
	return 0;
}

之后我们考虑化简一下式子

我们有fn-sigma(C(n,i)*fi*2)=0

之后构造多项式可以得到f - g*f =1 (因为f0=1)

则f = 1/(1-g)

多项式求逆即可

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cstdlib>
#include<algorithm>
#define G 3
using namespace std;

typedef long long LL;
const int mod=998244353;
const int maxn=300010;
int n,N,len,ans=0;
int jc[maxn],inv[maxn],rev[maxn];
int A[maxn],B[maxn];

int pow_mod(int v,int p){
	int tmp=1;
	while(p){
		if(p&1)tmp=1LL*tmp*v%mod;
		v=1LL*v*v%mod;p>>=1;
	}return tmp;
}
void FFT(int *A,int n,int flag){
	for(int i=0;i<n;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
	for(int k=2;k<=n;k<<=1){
		int mk=(k>>1);
		int wn=pow_mod(G,flag==1?(mod-1)/k:(mod-1)-(mod-1)/k);
		for(int i=0;i<n;i+=k){
			int w=1;
			for(int j=0;j<mk;++j){
				int a=i+j,b=i+j+mk;
				int x=A[a],y=1LL*A[b]*w%mod;
				A[a]=x+y;if(A[a]>=mod)A[a]-=mod;
				A[b]=x-y;if(A[b]<0)A[b]+=mod;
				w=1LL*w*wn%mod;
			}
		}
	}
	if(flag==-1){
		int c=pow_mod(n,mod-2);
		for(int i=0;i<n;++i)A[i]=1LL*A[i]*c%mod;
	}return;
}
void Get_inv(int n){
	if(n==1){
		B[0]=pow_mod(A[0],mod-2);
		return;
	}
	Get_inv(n>>1);
	int now=(n<<1);
	for(len=0;(1<<len)<now;++len);
	for(int i=0;i<now;++i)rev[i]=rev[i>>1]>>1|((i&1)<<(len-1));
	static int tmp[maxn];
	for(int i=0;i<n;++i)tmp[i]=A[i];
	FFT(tmp,now,1);FFT(B,now,1);
	for(int i=0;i<now;++i)B[i]=1LL*B[i]*(2LL-1LL*tmp[i]*B[i]%mod+mod)%mod;
	FFT(B,now,-1);
	memset(B+n,0,sizeof(int)*n);
}

int main(){
	freopen("heoi2016_sum.in","r",stdin);
	freopen("heoi2016_sum.out","w",stdout);
	scanf("%d",&n);
	for(N=1;N<=n;N<<=1);
	jc[0]=1;
	for(int i=1;i<N;++i)jc[i]=1LL*jc[i-1]*i%mod;
	inv[N-1]=pow_mod(jc[N-1],mod-2);
	for(int i=N-2;i>=0;--i)inv[i]=1LL*inv[i+1]*(i+1)%mod;
	for(int i=1;i<N;++i){
		A[i]=(inv[i]<<1);
		if(A[i]>=mod)A[i]-=mod;
	}A[0]--;if(A[0]<0)A[0]+=mod;
	Get_inv(N);
	for(int i=0;i<=n;++i){
		ans=ans+1LL*(mod-B[i])*jc[i]%mod;
		if(ans>=mod)ans-=mod;
	}printf("%d
",ans);
	return 0;
}

另外，求贝尔数的CDQ+FFT的程序

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#define G 11
using namespace std;

typedef long long LL;
const int maxn=300010;
const int mod=786433;
int n,N,len;
int jc[maxn],inv[maxn],rev[maxn];
int f[maxn],A[maxn],B[maxn],C[maxn];

int pow_mod(int v,int p){
	int tmp=1;
	while(p){
		if(p&1)tmp=1LL*tmp*v%mod;
		v=1LL*v*v%mod;p>>=1;
	}return tmp;
}
void FFT(int *A,int n,int flag){
	for(int i=0;i<n;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
	for(int k=2;k<=n;k<<=1){
		int mk=(k>>1);
		int wn=pow_mod(G,flag==1?(mod-1)/k:(mod-1)-(mod-1)/k);
		for(int i=0;i<n;i+=k){
			int w=1;
			for(int j=0;j<mk;++j){
				int a=i+j,b=i+j+mk;
				int x=A[a],y=1LL*A[b]*w%mod;
				A[a]=x+y;if(A[a]>=mod)A[a]-=mod;
				A[b]=x-y;if(A[b]<0)A[b]+=mod;
				w=1LL*w*wn%mod;
			}
		}
	}
	if(flag==-1){
		int c=pow_mod(n,mod-2);
		for(int i=0;i<n;++i)A[i]=1LL*A[i]*c%mod;
	}return;
}
void Solve(int L,int R){
	if(L==R)return;
	int mid=(L+R)>>1;
	Solve(L,mid);
	for(N=1,len=0;N<=(R-L+1);N<<=1,len++);
	A[0]=B[0]=0;
	for(int i=1;i<N;++i)A[i]=0,B[i]=inv[i-1];
	for(int i=L;i<=mid;++i)A[i-L]=1LL*f[i]*inv[i]%mod;
	if(R-L<=500){
		int lim=R-L+1;
		for(int i=0;i<lim;++i){
			C[i]=0;
			for(int j=0;j<=i;++j){
				C[i]=C[i]+1LL*A[j]*B[i-j]%mod;
				if(C[i]>=mod)C[i]-=mod;
			}
		}
	}else{
		for(int i=0;i<N;++i)rev[i]=rev[i>>1]>>1|((i&1)<<(len-1));
		FFT(A,N,1);FFT(B,N,1);
		for(int i=0;i<N;++i)C[i]=1LL*A[i]*B[i]%mod;
		FFT(C,N,-1);
	}
	for(int i=mid+1;i<=R;++i){
		f[i]=f[i]+1LL*jc[i-1]*C[i-L]%mod;
		if(f[i]>=mod)f[i]-=mod;
	}Solve(mid+1,R);
}

int main(){
	scanf("%d",&n);
	jc[0]=1;
	for(int i=1;i<=n;++i)jc[i]=1LL*jc[i-1]*i%mod;
	inv[n]=pow_mod(jc[n],mod-2);
	for(int i=n-1;i>=0;--i)inv[i]=1LL*inv[i+1]*(i+1)%mod;
	f[0]=1;Solve(0,n);
	printf("%d
",f[n]);
	return 0;
}

但是我并不知道怎么转成多项式求逆，如果有老司机知道还请告诉我QAQ

留下的一些坑：

1、城市规划的多项式求ln的写法

2、多项式开根

3、贝尔数的多项式求exp的写法

感觉很多CDQ+FFT都能很轻松的转成多项式求逆啊

虽然转了之后并没有快多少