多项式求逆 学习总结
Posted _Vertical
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了多项式求逆 学习总结相关的知识,希望对你有一定的参考价值。
感觉蒟蒻现在学多项式求逆貌似有些晚了
不过真的很有意思了(然而省选的时候自己还在玩泥巴什么也不会
多项式求逆是基于倍增的
假设我们知道
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\n",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\n",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\n",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\n",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\n",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\n",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\n",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\n",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\n",f[n]); return 0; }
但是我并不知道怎么转成多项式求逆,如果有老司机知道还请告诉我QAQ
留下的一些坑:
1、城市规划的多项式求ln的写法
2、多项式开根
3、贝尔数的多项式求exp的写法
感觉很多CDQ+FFT都能很轻松的转成多项式求逆啊
虽然转了之后并没有快多少
以上是关于多项式求逆 学习总结的主要内容,如果未能解决你的问题,请参考以下文章