BZOJ 3684 大朋友和多叉树 FFT+拉格朗日反演

题目大意：给定$n$和集合$S$，求满足下列要求的多叉树的个数：
1.每个非叶节点的子节点数量在集合$S$中
2.每个叶节点的权值为$1$，每个非叶节点的权值为子节点权值之和
3.根节点的权值为$n$
注意每个节点的子节点有顺序

令$f_i$表示根节点权值为$i$的神犇二叉树个数，$F(x)$为$f_i$的生成函数，$C(x)$为$S$的生成函数，那么有：
$F(x)=\sum_{i\in S}F^i(x)+x$
$F(x)=C(F(x))+x$
$F(x)-C(F(x))=x$
不妨令$G(x)=1-C(x)$
那么有：
$G(F(x))=x$
因此$F(x)$是$G(x)$的复合逆
拉格朗日反演：
$[x^n]F(x)=\frac1n[x^{n-1}](\frac x{G(x)})^n$
其中$[x^n]F(x)$表示$F(x)$的$n$次项系数
(我并不知道怎么证明，去问策爷吧
然后就是FFT一通乱求的事了……
后面那个$n$次方用$ln+exp$实测比快速幂快
然后T了……
孩子FFT老是T，多半是废了！
于是怒换wyfcyx的非递归FFT模板
结果换模板之前极限数据16s，换完后极限数据5s
[摔西瓜]坑爹呢么这不是！
我说我怎么FFT写一道T一道！
妈妈再也不用担心我的FFT了！
关于此题的其他细节敬请参阅金策在Wc2015上的营员交流《生成函数的运算与应用》

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define M 263000
#define MOD 950009857
#define G 7
using namespace std;
int n,m,d;
long long inv[M];
void Linear_Shaker()
{
    int i;
    for(inv[1]=1,i=2;i<=d<<1;i++)
        inv[i]=(MOD-MOD/i)*inv[MOD%i]%MOD;
}
long long Quick_Power(long long x,int y)
{
    long long re=1;
    while(y)
    {
        if(y&1) (re*=x)%=MOD;
        (x*=x)%=MOD; y>>=1;
    }
    return re;
}
/*
void FFT(int a[],int n,int type)
{
    static int temp[M];
    int i;
    if(n==1) return ;
    for(i=0;i<n;i+=2)
        temp[i>>1]=a[i],temp[i+n>>1]=a[i+1];
    memcpy(a,temp,sizeof(a[0])*n);
    int *l=a,*r=a+(n>>1);
    FFT(l,n>>1,type);
    FFT(r,n>>1,type);
    long long w=Quick_Power(G,(long long)(MOD-1)/n*type%(MOD-1)),wn=1;
    for(i=0;i<n>>1;i++,(wn*=w)%=MOD)
        temp[i]=(l[i]+wn*r[i])%MOD,temp[i+(n>>1)]=(l[i]-wn*r[i]%MOD+MOD)%MOD;
    memcpy(a,temp,sizeof(a[0])*n);
}
*/
inline int Revbit(int x,int bit){
    int r=0;
    for(int i=0;i<bit;++i)
        if((x>>i)&1)
            r+=1<<(bit-1-i);
    return r;
}
inline void FFT(int a[],int n,int rev){
    int i,j,k,w,wn,t,bit=0;
    for(int tmp=n;tmp^1;tmp>>=1) ++bit;
    static int b[M];
    for(i=0;i<n;i++)
        b[Revbit(i,bit)]=a[i];
    memcpy(a,b,sizeof(a[0])*n);
    for(k=2;k<=n;k<<=1)
        for(wn=Quick_Power(G,(long long)(MOD-1)/k*rev%(MOD-1)),i=0;i<n;i+=k)
            for(w=1,j=0;j<k/2;++j,w=(long long)w*wn%MOD)
                t=(long long)w*a[i+j+k/2]%MOD,a[i+j+k/2]=(a[i+j]+MOD-t)%MOD,a[i+j]=(a[i+j]+t)%MOD;  
}
void Get_Inv(int a[],int b[],int n)
{
    static int temp[M];
    int i;
    if(n==1)
    {
        b[0]=Quick_Power(a[0],MOD-2);
        return ;
    }
    Get_Inv(a,b,n>>1);
    memcpy(temp,a,sizeof(a[0])*n);
    memset(temp+n,0,sizeof(a[0])*n);
    FFT(temp,n<<1,1);
    FFT(b,n<<1,1);
    for(i=0;i<n<<1;i++)
        temp[i]=(long long)b[i]*(2-(long long)temp[i]*b[i]%MOD+MOD)%MOD;
    FFT(temp,n<<1,MOD-2);
    for(i=0;i<n;i++)
        b[i]=(long long)temp[i]*inv[n<<1]%MOD;
    memset(b+n,0,sizeof(a[0])*n);
}
/*
void Quick_Power(int x[],int re[],int y)
{
    int i;
    long long inv=Quick_Power(d<<1,MOD-2);
    re[0]=1;
    while(y)
    {
        FFT(x,d<<1,1);

        if(y&1)
        {
            FFT(re,d<<1,1);
            for(i=0;i<d<<1;i++)
                re[i]=(long long)x[i]*re[i]%MOD;
            FFT(re,d<<1,MOD-2);
            for(i=0;i<d;i++)
                re[i]=inv*re[i]%MOD;
            memset(re+d,0,sizeof(re[0])*d);
        }

        for(i=0;i<d<<1;i++)
            x[i]=(long long)x[i]*x[i]%MOD;
        FFT(x,d<<1,MOD-2);
        for(i=0;i<d;i++)
            x[i]=inv*x[i]%MOD;
        memset(x+d,0,sizeof(x[0])*d);

        y>>=1;
    }
}
*/
void Get_Ln(int a[],int b[],int n)
{
    static int a_[M],a_inv[M];
    int i;
    Get_Inv(a,a_inv,n);
    for(i=0;i<n-1;i++)
        a_[i]=(long long)a[i+1]*(i+1)%MOD;
    FFT(a_,n<<1,1);
    FFT(a_inv,n<<1,1);
    for(i=0;i<n<<1;i++)
        b[i]=(long long)a_[i]*a_inv[i]%MOD*inv[n<<1]%MOD;
    FFT(b,n<<1,MOD-2);
    for(i=n-1;i;i--)
        b[i]=b[i-1]*inv[i]%MOD;
    b[0]=0;
    memset(b+n,0,sizeof(b[0])*n);
    memset(a_,0,sizeof(a_[0])*n<<1);
    memset(a_inv,0,sizeof(a_inv[0])*n<<1);
}
void Get_Exp(int a[],int b[],int n)
{
    static int temp[M];
    int i;
    if(n==1)
    {
        b[0]=1;
        return ;
    }
    Get_Exp(a,b,n>>1);
    memset(temp,0,sizeof(temp[0])*n<<1);
    Get_Ln(b,temp,n);
    for(i=0;i<n;i++)
        temp[i]=((i==0)+MOD-temp[i]+a[i])%MOD;
    FFT(temp,n<<1,1);
    FFT(b,n<<1,1);
    for(i=0;i<n<<1;i++)
        b[i]=(long long)b[i]*temp[i]%MOD;
    FFT(b,n<<1,MOD-2);
    for(i=0;i<n;i++)
        b[i]=inv[n<<1]*b[i]%MOD;
    memset(b+n,0,sizeof(b[0])*n);
}
int main()
{
    static int A[M],G_[M],G_ln[M],G_n[M];
    int i,x;
    cin>>n>>m;
    for(d=1;d<=n;d<<=1);
    Linear_Shaker();
    for(A[0]=1,i=1;i<=m;i++)
    {
        scanf("%d",&x);
        A[x-1]=MOD-1;
    }
    Get_Inv(A,G_,d);
    //Quick_Power(_G,G_d,n);
    Get_Ln(G_,G_ln,d);
    for(i=0;i<d;i++)
        G_ln[i]=(long long)G_ln[i]*n%MOD;
    Get_Exp(G_ln,G_n,d);

    //static int test[M];
    //Get_Ln(G_n,test,d);

    cout<<G_n[n-1]*Quick_Power(n,MOD-2)%MOD<<endl;
    return 0;
}