多项式求逆 学习总结
感觉蒟蒻现在学多项式求逆貌似有些晚了
不过真的很有意思了(然而省选的时候自己还在玩泥巴什么也不会
多项式求逆是基于倍增的
假设我们知道
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都能很轻松的转成多项式求逆啊
虽然转了之后并没有快多少
