BZOJ 3684 大朋友和多叉树（生成函数+FFT）

Description

我们的大朋友很喜欢计算机科学，而且尤其喜欢多叉树。对于一棵带有正整数点权的有根多叉树，如果它满足这样的性质，我们的大朋友就会将其称作神犇的：点权为1的结点是叶子结点；对于任一点权大于1的结点u，u的孩子数目deg[u]属于集合D，且u的点权等于这些孩子结点的点权之和。

给出一个整数s，你能求出根节点权值为s的神犇多叉树的个数吗？请参照样例以更好的理解什么样的两棵多叉树会被视为不同的。

我们只需要知道答案关于950009857(453∗221+1，一个质数）取模后的值。

Input

第一行有2个整数s,m。

第二行有m个互异的整数，d[1],d[2],…,d[m]，为集合D中的元素。

Output

输出一行仅一个整数，表示答案模950009857的值。

Sample Input

4 2

2 3

Sample Output

10

Solution

设满足条件的树的生成函数为F(x)，则有F(x)=x+∑c∈D(F(x))c

设A(x)=x−∑c∈Dxc，则A(F(x))=x，由拉格朗日反演，[xn]F(x)=1n[xn−1](xA(x))n

而(xA(x))n=e−nln(A(x)x)，故只需多项式A(x)x取对数后乘上−n然后取指数即可

Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
#define maxn 100005
#define maxfft 262144+5
#define mod 950009857
const double pi=acos(-1.0);
struct cp
{
double a,b;
cp operator +(const cp &o)const {return (cp){a+o.a,b+o.b};}
cp operator -(const cp &o)const {return (cp){a-o.a,b-o.b};}
cp operator *(const cp &o)const {return (cp){a*o.a-b*o.b,b*o.a+a*o.b};}
cp operator *(const double &o)const {return (cp){a*o,b*o};}
cp operator !() const{return (cp){a,-b};}
}w[maxfft];
int pos[maxfft];
void fft_init(int len)
{
int j=0;
while((1<<j)<len)j++;
j--;
for(int i=0;i<len;i++)
pos[i]=pos[i>>1]>>1|((i&1)<<j);
}
void fft(cp *x,int len,int sta)
{
for(int i=0;i<len;i++)
if(i<pos[i])swap(x[i],x[pos[i]]);
w[0]=(cp){1,0};
for(unsigned i=2;i<=len;i<<=1)
{
cp g=(cp){cos(2*pi/i),sin(2*pi/i)*sta};
for(int j=i>>1;j>=0;j-=2)w[j]=w[j>>1];
for(int j=1;j<i>>1;j+=2)w[j]=w[j-1]*g;
for(int j=0;j<len;j+=i)
{
cp *a=x+j,*b=a+(i>>1);
for(int l=0;l<i>>1;l++)
{
cp o=b[l]*w[l];
b[l]=a[l]-o;
a[l]=a[l]+o;
}
}
}
if(sta==-1)for(int i=0;i<len;i++)x[i].a/=len,x[i].b/=len;
}
cp x[maxfft],y[maxfft],z[maxfft];
int temp[maxfft];
void FFT(int *a,int *b,int n,int m,int *c)
{
if(n<=100&&m<=100||min(n,m)<=5)
{
for(int i=0;i<n+m-1;i++)temp[i]=0;
for(int i=0;i<n;i++)
for(int j=0;j<m;j++)
{
temp[i+j]+=(ll)a[i]*b[j]%mod;
if(temp[i+j]>=mod)temp[i+j]-=mod;
}
for(int i=0;i<n+m-1;i++)c[i]=temp[i];
return ;
}
int len=1;
while(len<n+m)len<<=1;
fft_init(len);
for(int i=0;i<len;i++)
{
int aa=i<n?a[i]:0,bb=i<m?b[i]:0;
x[i]=(cp){(aa>>15),(aa&32767)},y[i]=(cp){(bb>>15),(bb&32767)};
}
fft(x,len,1),fft(y,len,1);
for(int i=0;i<len;i++)
{
int j=len-1&len-i;
z[i]=((x[i]+!x[j])*(y[i]-!y[j])+(x[i]-!x[j])*(y[i]+!y[j]))*(cp){0,-0.25};
}
fft(z,len,-1);
for(int i=0;i<n+m-1;i++)
{
ll ta=(ll)(z[i].a+0.5)%mod;
ta=(ta<<15)%mod;
c[i]=ta;
}
for(int i=0;i<len;i++)
{
int j=len-1&len-i;
z[i]=(x[i]-!x[j])*(y[i]-!y[j])*(cp){-0.25,0}+(x[i]+!x[j])*(y[i]+!y[j])*(cp){0,0.25};
}
fft(z,len,-1);
for(int i=0;i<n+m-1;i++)
{
ll ta=(ll)(z[i].a+0.5)%mod,tb=(ll)(z[i].b+0.5)%mod;
ta=(ta+(tb<<30))%mod;
c[i]=(c[i]+ta)%mod;
}
}
int inv[maxn];
void init(int n=100001)
{
inv[1]=1;
for(int i=2;i<=n;i++)inv[i]=mod-(ll)(mod/i)*inv[mod%i]%mod;
}
int temp1[maxfft],temp2[maxfft],temp3[maxfft],temp4[maxfft];
void Poly_Inv(int *poly,int n,int *ans)
{
ans[0]=inv[poly[0]];
for(int i=2;i<=n;i<<=1)
{
FFT(poly,ans,i,i/2,temp1);
FFT(ans,temp1+i/2,i/2,i/2,temp1);
for(int j=0;j<i/2;j++)ans[j+i/2]=temp1[j]==0?0:mod-temp1[j];
}
}
void Poly_Log(int *poly,int n,int *ans)
{
Poly_Inv(poly,n,temp2);
for(int i=0;i<n-1;i++)ans[i]=(ll)poly[i+1]*(i+1)%mod;
FFT(ans,temp2,n-1,n,ans);
for(int i=n-1;i>0;i--)ans[i]=(ll)ans[i-1]*inv[i]%mod;
ans[0]=0;
}
void Poly_Exp(int *poly,int n,int *ans)
{
if(n==1)
{
ans[0]=1;
return ;
}
Poly_Exp(poly,n/2,ans);
Poly_Log(ans,n,temp3);
for(int i=0;i<n;i++)
{
temp3[i]=poly[i]-temp3[i];
if(temp3[i]<0)temp3[i]+=mod;
}
temp3[0]++;
if(temp3[0]==mod)temp3[0]=0;
FFT(ans,temp3,n,n,ans);
for(int i=n;i<2*n;i++)ans[i]=0;
}
void Poly_Root(int *poly,int n,int *ans)
{
ans[0]=1;
for(int i=2;i<=n;i<<=1)
{
Poly_Inv(ans,i,temp4);
FFT(ans,ans,i/2,i/2,ans);
for(int j=0;j<i;j++)ans[j]=(ll)(ans[j]+poly[j])*inv[2]%mod;
FFT(ans,temp4,i,i,ans);
for(int j=i;j<2*i;j++)ans[j]=0;
}
}
int s,m,f[maxfft],g[maxfft];
int main()
{
init();
while(~scanf("%d%d",&s,&m))
{
memset(f,0,sizeof(f));
while(m--)
{
int temp;
scanf("%d",&temp);
f[temp-1]=mod-1;
}
f[0]=1;
int len=1;
while(len<=s)len<<=1;
Poly_Log(f,len,g);
for(int i=0;i<len;i++)g[i]=(mod-(ll)s*g[i]%mod)%mod;
Poly_Exp(g,len,f);
printf("%d\n",(ll)inv[s]*f[s-1]%mod);
}
return 0;
}