大朋友与多叉树

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

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

我们只需要知道答案关于 950009857（\(453 \times 2^{21} + 1\)，一个质数）取模后的值。

题解

https://jkloverdcoi.github.io/2019/12/01/%E6%8B%89%E6%A0%BC%E6%9C%97%E6%97%A5%E5%8F%8D%E6%BC%94%E5%AD%A6%E4%B9%A0%E7%AC%94%E8%AE%B0/
https://www.cnblogs.com/xiao-ju-ruo-xjr/p/8455362.html

我们还从暴力DP入手。令 f(i) 表示根节点权值为 i 的树有多少个，令 g(i, j) 表示由 j 棵树组成的根节点权值总和为 i 的森林有多少种，则显然由以下转移：

\[ f(i) = \sum_{j \in D} g(i, j) \\ f(1) = 1 \\ g(i, j) = \sum_{k = 1}^j g(i - k, j - 1) \cdot f(k) \]

现在如果用 F(x) 表示 f(i) 的生成函数，\(G_j(x)\) 表示 g(i, j) 的生成函数，不难发现 \(G_j(x) = F(x)^j\)，那么就能得到这样一个方程：

\[ F(x) = x + \sum_{j \in D} F(x)^j \]

为什么要加一个 x 呢？从DP的边界条件看来 f(1) = 1 而 f(0) = 0；此外根据题意 D 中元素都 \(\ge 2\)，而显然 F(x) 的常数项为 0，那么 \(\sum_{j \in D} F(x)^j\) 的最低次项就是 D 中最小元素了，显然比一次项要高，所以我们需要加一个 x。

但是这个方程我并不会解，怎么办？这就是拉格朗日反演定理的用处了。先看一下这个定理在什么条件下能用。

f(x) 与 g(x) 常数项为 0，且 g(f(x)) = x，那么若已知 g(x)，我们可以快速求 f(x) 的某一项（不妨设要求的是 \(x^n\) 项），我们有：

\[ [x^n]f(x) = \frac{1}{n} [x^{-1}] \left( \frac{1}{g(x)} \right) ^n \]

事实上若满足常数项都是 0 的条件，f(g(x)) = g(f(x)) = x 是成立的，所以网上对这个定理的描述可能稍有不同，其实都是等价的。

可能对上面的公式你还不明白 \([x^{-1}]\) 是什么操作。注意到 g(x) 的常数项为 0，所以它的逆元是不存在的，即没有任何整式能够表示 \(\frac{1}{g(x)}\)，这时候就需要扩充一下这个域，引入分式域。在这里所有的多项式能够被表示成 \(\cdots + a_{-2} x^{-2} + a_{-1} x^{-1} + a_0 + a_1 x + a_2 x^2 + \cdots\)，可以证明这种形式能够表示所有的 \(\frac{A(x)}{B(x)}\)，其中 A(x) 和 B(x) 都是整式。说白了 \(\frac{1}{g(x)}\) 是一个分式，所以它有 \(x^{-1}\) 项。

但是你还想知道如何求一个分式的 \(x^{-1}\) 项。其实我们不需要真正求一个分式，我们可以把上面的 g(x) 的末尾的系数为 0 的项都去掉，然后正常地求逆元，然后再左移对应项数。即将上式变为：

\[ [x^n]f(x) = \frac{1}{n} [x^{dn-1}] \left( \frac{x^d}{g(x)} \right) ^n \]

其中 \(d\) 表示 \(g(x)\) 的系数中前缀 \(0\) 的个数。这个式子就是求一下 \(\mod x^{dn}\) 下 g'(x) 的逆元就好了，g'(x) 是 g(x) 左移直到常数项非零后的多项式。

回到这题，不难发现 (1) 式移项即可得到：

\[ F(x) - \sum_{j \in D} F(x)^j = x \\ G(F(x)) = x \]

所以得到了 \(G(x) = x - \sum_{j \in D} x^j\) 这个多项式，问题就解决了。

时间复杂度 \(O(n \log n)\)。

CO int N=262144;
int omg[2][N],rev[N];
int fac[N],inv[N],ifac[N];

void NTT(poly&a,int dir){
    int lim=a.size(),len=log2(lim);
    for(int i=0;i<lim;++i) rev[i]=rev[i>>1]>>1|(i&1)<<(len-1);
    for(int i=0;i<lim;++i)if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int i=1;i<lim;i<<=1)
        for(int j=0;j<lim;j+=i<<1)for(int k=0;k<i;++k){
            int t=mul(omg[dir][N/(i<<1)*k],a[j+i+k]);
            a[j+i+k]=add(a[j+k],mod-t),a[j+k]=add(a[j+k],t);
        }
    if(dir==1){
        int ilim=fpow(lim,mod-2);
        for(int i=0;i<lim;++i) a[i]=mul(a[i],ilim);
    }
}
poly operator~(poly a){
    int n=a.size();
    poly b={fpow(a[0],mod-2)};
    if(n==1) return b;
    a.resize(1<<(int)ceil(log2(n)));
    for(int lim=2;lim<2*n;lim<<=1){
        poly c(a.begin(),a.begin()+lim);
        c.resize(lim<<1),NTT(c,0);
        b.resize(lim<<1),NTT(b,0);
        for(int i=0;i<lim<<1;++i) b[i]=mul(2+mod-mul(c[i],b[i]),b[i]);
        NTT(b,1),b.resize(lim);
    }
    return b.resize(n),b;
}
poly log(poly a){
    int n=a.size();
    poly b=~a;
    for(int i=0;i<n-1;++i) a[i]=mul(a[i+1],i+1);
    a.resize(n-1);
    int lim=1<<(int)ceil(log2(2*n-2));
    a.resize(lim),NTT(a,0);
    b.resize(lim),NTT(b,0);
    for(int i=0;i<lim;++i) a[i]=mul(a[i],b[i]);
    NTT(a,1),a.resize(n);
    for(int i=n-1;i>=1;--i) a[i]=mul(a[i-1],inv[i]);
    return a[0]=0,a;
}
poly exp(poly a){
    int n=a.size();
    poly b={1}; // a[0]=0
    if(n==1) return b;
    a.resize(1<<(int)ceil(log2(n)));
    for(int lim=2;lim<2*n;lim<<=1){
        b.resize(lim);poly c=log(b);
        c[0]=add(1+a[0],mod-c[0]);
        for(int i=1;i<lim;++i) c[i]=add(a[i],mod-c[i]);
        c.resize(lim<<1),NTT(c,0);
        b.resize(lim<<1),NTT(b,0);
        for(int i=0;i<lim<<1;++i) b[i]=mul(c[i],b[i]);
        NTT(b,1),b.resize(lim);
    }
    return b.resize(n),b;
}

int main(){
    omg[0][0]=1,omg[0][1]=fpow(7,(mod-1)/N);
    omg[1][0]=1,omg[1][1]=fpow(omg[0][1],mod-2);
    fac[0]=fac[1]=1;
    inv[0]=inv[1]=1;
    ifac[0]=ifac[1]=1;
    for(int i=2;i<N;++i){
        omg[0][i]=mul(omg[0][i-1],omg[0][1]);
        omg[1][i]=mul(omg[1][i-1],omg[1][1]);
        fac[i]=mul(fac[i-1],i);
        inv[i]=mul(mod-mod/i,inv[mod%i]);
        ifac[i]=mul(ifac[i-1],inv[i]);
    }
    int n=read<int>();
    poly g(n);
    for(int m=read<int>();m--;){
        int d=read<int>();
        if(d-1<n) g[d-1]=mod-1;
    }
    g[0]=1;
    g=log(~g);
    for(int i=0;i<n;++i) g[i]=mul(g[i],n);
    g=exp(g);
    int ans=mul(g[n-1],inv[n]);
    printf("%d\n",ans);
    return 0;
}