「bzoj3684」 大朋友与多叉树

当年这题和黄队一起看过，但是没做，现在只剩我自己了（结果黄队一人独自把这题 A 了

Description

我们的大朋友很喜欢计算机科学，而且尤其喜欢多叉树。对于一棵带有正整数点权的有根多叉树，如果它满足这样的性质，我们的大朋友就会将其称作神犇的：点权为 $1$ 的结点是叶子结点；对于任一点权大于 $1$ 的结点 $u$，$u$ 的孩子数目$deg[u]$ 属于集合 $D$ ，且 $u$ 的点权等于这些孩子结点的点权之和。
给出一个整数 $s$，你能求出根节点权值为 $s$ 的神犇多叉树的个数吗？请参照样例以更好的理解什么样的两棵多叉树会被视为不同的。
我们只需要知道答案关于 $950009857$（ $453\times2^{21}+1$，一个质数）取模后的值。

Solution

这个题也是个裸题8，但是是第一次做，于是就写了一篇博客。

首先考虑一个暴力的 DP ：$f_i$ 表示叶子节点数为 $i$ 的合法多叉树数量，转移时枚举自己有几个儿子，做一个类似背包的转移 ：$f_i = \sum_{j\in D}\left((\prod_{k=1}^jf_{son_i})\times[\sum_{k=1}^json_i = i]\right)$，因为一个点自己就是叶子，所以 $f_1 = 1$。

我们令 $F(x)$ 为 $f$ 的生成函数，根据上面那个式子，我们可以写出 $F(x) = \sum_{i\in D}F^i(x) + x$ ，简单移一下项，得到 $F(x) - \sum_{i\in D} F^i(x) = x$ ，毛估估一下，$F(x)$ 的复合逆 $G(x)$ 就是 $x - \sum_{i\in D} x^i$，题目要我们求的东西是 $[x^s]F(x)$，套一个拉格朗日反演上去，$[x^n]F(x) = \frac1n[x^{-1}]\left( \frac{1}{G(x)} \right)^n$ ，大概就做完了。

code

#include <map>
#include <set>
#include <ctime>
#include <queue>
#include <stack>
#include <cmath>
#include <vector>
#include <bitset>
#include <cstdio>
#include <cctype>
#include <string>
#include <numeric>
#include <cstring>
#include <cassert>
#include <climits>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#include <functional>
using namespace std ;
#define rep(i, a, b) for (register int i = (a); i <= (b); ++i)
#define per(i, a, b) for (register int i = (a); i >= (b); --i)
#define loop(it, v) for (auto it = v.begin(); it != v.end(); it++)
#define cont(i, x) for (register int i = head[x]; i; i = edge[i].nex)
#define clr(a) memset(a, 0, sizeof(a))
#define ass(a, cnt) memset(a, cnt, sizeof(a))
#define cop(a, b) memcpy(a, b, sizeof(a))
#define lowbit(x) (x & -x)
#define all(x) x.begin(), x.end()
#define SC(t, x) static_cast <t> (x)
#define ub upper_bound
#define lb lower_bound
#define pqueue priority_queue
#define mp make_pair
#define pb push_back
#define pof pop_front
#define pob pop_back
#define fi first
#define se second
#define y1 y1_
#define Pi acos(-1.0)
#define iv inline void
#define enter putchar('\n')
#define siz(x) ((int)x.size())
#define file(x) freopen(x".in", "r", stdin),freopen(x".out", "w", stdout)
typedef double db ;
typedef long long ll ;
typedef unsigned long long ull ;
typedef pair <int, int> pii ;
typedef vector <int> vi ;
typedef vector <pii> vii ;
typedef queue <int> qi ;
typedef queue <pii> qii ;
typedef set <int> si ;
typedef map <int, int> mii ;
typedef map <string, int> msi ;
const int maxn = 1e6 + 100 ;
const int inf = 0x3f3f3f3f ;
const int iinf = 1 << 30 ;
const ll linf = 2e18 ;
const int mod = 950009857 ;
const double eps = 1e-5 ;
template <class T = ll> T read()
{
    T f = 1, a = 0;
    char ch = getchar() ;
    while (!isdigit(ch)) { if (ch == '-') f = -1 ; ch = getchar() ; }
    while (isdigit(ch)) { a =  ((a << 3) + (a << 1) + ch - '0') % mod; ch = getchar() ; }
    return a * f ;
}

namespace polynomial
{
    inline int power(int a, int b)
	{
        register int ret = 1;
        for(; b; b >>= 1, a = 1ll * a * a % mod) if(b & 1) ret = 1ll * ret * a % mod;
        return ret;
    }

    int rev[maxn], g = 7, gg[2][maxn], lim_inv[maxn];

    inline void init()
	{
        for(register int j = 1; j < (maxn >> 1); j <<= 1)
        {
            register int g1 = power(g, (mod - 1) / j), g0 = power(g1, mod - 2);
            lim_inv[j] = power(j, mod - 2);
            gg[0][j] = 1, gg[1][j] = 1;
            rep(i, 1, j - 1) gg[0][j + i] = 1ll * gg[0][j + i - 1] * g0 % mod, gg[1][j + i] = 1ll * gg[1][j + i - 1] * g1 % mod;
        }
    }
    
    inline void ntt(int *a, int lim, int opt)
	{
        rep(i, 0, lim - 1) if(i < rev[i]) swap(a[i], a[rev[i]]);
        for(register int i = 1; i < lim; i <<= 1)
            for(register int j = 0; j < lim; j += (i << 1))
                for(register int k = 0; k < i; ++ k)
                {
                    register int x = a[j + k], y = 1ll * a[i + j + k] * gg[opt][(i << 1) + k] % mod;
                    a[j + k] = (x + y) % mod;
                    a[i + j + k] = (x + mod - y) % mod;
                }
        if(!opt) rep(i, 0, lim - 1) a[i] = 1ll * a[i] * lim_inv[lim] % mod;
    }

    inline void mul(int *a, int *b, int len)
	{
        register int lim = 1;
        while(lim < (len << 1)) lim <<= 1;
        rep(i, 0, lim - 1) rev[i] = (rev[i >> 1] >> 1) + (i & 1) * (lim >> 1);  
        ntt(a, lim, 1), ntt(b, lim, 1);
        rep(i, 0, lim - 1) a[i] = 1ll * a[i] * b[i] % mod;
        ntt(a, lim, 0);
    }

    inline void poly_mul(int *a, int *b, int len)
	{
        register int lim = 1;
        while(lim < (len << 1)) lim <<= 1;
        rep(i, 0, lim - 1) rev[i] = (rev[i >> 1] >> 1) + (i & 1) * (lim >> 1);
        // puts("poly : mul");
        // rep(i, 0, lim - 1) printf("%d %d\n", a[i], b[i]);
        ntt(a, lim, 1), ntt(b, lim, 1);
        rep(i, 0, lim - 1) a[i] = 1ll * a[i] * b[i] % mod;
        ntt(a, lim, 0);
        rep(i, len, lim) a[i] = b[i] = 0;
    }

    int ta[maxn], tb[maxn];

    void dac_mul(int l, int r, int *f, int *g)
	{
        if(l == r) return ;
        int mid = (l + r) >> 1, len = r - l + 1, lim = 1;
        dac_mul(l, mid, f, g);
        while(lim < (len << 1)) lim <<= 1;
        rep(i, 0, lim - 1) ta[i] = tb[i] = 0;
        rep(i, 0, mid - l) ta[i] = f[i + l];
        rep(i, 0, r - l - 1) tb[i] = g[i + 1];
        mul(ta, tb, len);
        rep(i, mid + 1, r) (f[i] += ta[i - l - 1]) %= mod;
        dac_mul(mid + 1, r, f, g);
    }

    inline int num_plus(int &a, int b)
	{
        return (a += b) >= mod ? a -= mod : a;
    }

    inline void minus(int *f, int *g, int len)
	{
        rep(i, 0, len - 1) num_plus(f[i], mod - g[i]);
        return ;
    }

    inline void plus(int *f, int *g, int len)
	{
        rep(i, 0, len - 1) num_plus(f[i], g[i]);
        return ;
    }

    inline void num_mul(int *f, int x, int len)
	{
        rep(i, 0, len - 1) f[i] = 1ll * f[i] * x % mod;
        return ;
    }

    inline void intergral(int *f, int len)
	{
        per(i, len, 1) f[i] = 1ll * f[i - 1] * power(i, mod - 2) % mod;
        f[0] = 0;
    }

    inline void derivative(int *f, int len)
	{
        rep(i, 0, len - 1) f[i] = 1ll * f[i + 1] * (i + 1) % mod;
        return ;
    }

    int inv_ta[maxn], inv_tb[maxn];

    void inversion(int *f, int *a, int len)
	{
        rep(i, 0, len << 1) f[i] = 0;
        if(len == 1) return (void)(f[0] = power(a[0], mod - 2));
        inversion(f, a, (len + 1) >> 1);
        rep(i, 0, len - 1) inv_ta[i] = inv_tb[i] = f[i], (f[i] *= 2) %= mod;
        poly_mul(inv_ta, inv_tb, len);
        rep(i, 0, len - 1) inv_tb[i] = a[i];
        poly_mul(inv_ta, inv_tb, len), minus(f, inv_ta, len);
        rep(i, 0, len - 1) inv_ta[i] = inv_tb[i] = 0;
    }

    int sqrt_ta[maxn], sqrt_tb[maxn], sqrt_tc[maxn];

    void poly_sqrt(int *f, int *a, int len)
	{
        rep(i, 0, len << 1) f[i] = 0;
        if(len == 1) return (void)(f[0] = 1);
        poly_sqrt(f, a, (len + 1) >> 1);
        rep(i, 0, len - 1) sqrt_ta[i] = sqrt_tb[i] = f[i];
        poly_mul(sqrt_ta, sqrt_tb, len), plus(sqrt_ta, a, len);
        rep(i, 0, len - 1) sqrt_tb[i] = f[i] * 2 % mod, sqrt_tc[i] = 0;
        inversion(sqrt_tc, sqrt_tb, len), poly_mul(sqrt_ta, sqrt_tc, len);
        rep(i, 0, len - 1) f[i] = sqrt_ta[i], sqrt_ta[i] = sqrt_tb[i] = sqrt_tc[i] = 0;
    }

    int der[maxn], ln_inv[maxn];

    void poly_ln(int *f, int len)
	{
        rep(i, 0, len - 1) der[i] = ln_inv[i] = f[i];
        derivative(der, len), inversion(f, ln_inv, len), poly_mul(f, der, len), intergral(f, len);
        rep(i, 0, len - 1) der[i] = ln_inv[i] = 0;

    }

    int exp_ln[maxn], exp_mul[maxn];

    void poly_exp(int *f, int *a, int len)
	{
        rep(i, 0, len << 1) f[i] = 0;
        if(len == 1) return (void)(f[0] = 1);
        poly_exp(f, a, (len + 1) >> 1);
        rep(i, 0, len - 1) exp_ln[i] = f[i];
        poly_ln(exp_ln, len);
        exp_mul[0] = 1, minus(exp_mul, exp_ln, len), plus(exp_mul, a, len), poly_mul(f, exp_mul, len);
        rep(i, 0, len - 1) exp_ln[i] = exp_mul[i] = 0;
    }
    
    int prod[maxn];

    void poly_pow(int *f, int b, int len)
	{
        poly_ln(f, len), num_mul(f, b, len);
        rep(i, 0, len - 1) prod[i] = f[i];
        poly_exp(f, prod, len);
        rep(i, 0, len - 1) prod[i] = 0;
    }

    inline void poly_reverse(int *a, int len)
	{
        reverse(a, a + len);
    }

    int div_inv[maxn];

    inline void poly_div(int *f, int *a, int *b, int len_a, int len_b)
	{
        register int len = len_a - len_b + 1;
        if(len < 1) return ;
        poly_reverse(a, len_a), poly_reverse(b, len_b), inversion(div_inv, b, len);
        rep(i, 0, len - 1) f[i] = a[i];
        poly_mul(f, div_inv, len), poly_reverse(f, len), poly_reverse(a, len_a), poly_reverse(b, len_b);
        rep(i, 0, len - 1) div_inv[i] = 0;
    }

    int mod_div[maxn];

    inline void poly_mod(int *f, int *a, int *b, int len_a, int len_b)
	{
        rep(i, 0, len_b - 1) f[i] = b[i];
        poly_div(mod_div, a, b, len_a, len_b);
        poly_mul(f, mod_div, len_a), num_mul(f, mod - 1, len_a), plus(f, a, len_a);
        rep(i, 0, len_a) mod_div[i] = 0;
    }

    int mpe_mul[maxn * 10], mpe_begin[maxn], cnt, mpe_tmp[maxn];

    void poly_dac_mul(int *a, int l, int r, int tot)
	{
        if(l == r)
        {
            mpe_begin[tot] = ++ cnt;
            mpe_mul[cnt] = mod - a[l], mpe_mul[++ cnt] = 1;
            return ;
        }
        register int mid = (l + r) >> 1, len = r - l + 2;
        poly_dac_mul(a, l, mid, tot << 1), poly_dac_mul(a, mid + 1, r, tot << 1 | 1);
        mpe_mul[mpe_begin[tot] = ++ cnt] = 0;
        rep(i, 1, len - 1) mpe_mul[++ cnt] = 0;
        rep(i, 0, mid - l + 1) mpe_mul[mpe_begin[tot] + i] = mpe_mul[mpe_begin[tot << 1] + i];
        rep(i, 0, r - mid) mpe_tmp[i] = mpe_mul[mpe_begin[tot << 1 | 1] + i];
        poly_mul(mpe_mul + mpe_begin[tot], mpe_tmp, len);
        rep(i, 0, len - 1) mpe_tmp[i] = 0;
    }

    inline int poly_spe(int *f, int len, int point)
	{
        register int ret = 0, base = 1;
        rep(i, 0, len - 1) num_plus(ret, 1ll * base * f[i] % mod), base = 1ll * base * point % mod;
        return ret;
    }

    int mpe_r[maxn * 10], mpe_r_begin[maxn];

    const int shit = 256;

    void poly_mpe_done(int *f, int *a, int *b, int len, int l, int r, int tot)
	{
        register int last = len;
        len = r - l + 1;
        mpe_r[mpe_r_begin[tot] = ++ cnt] = 0;
        rep(i, 1, len - 1) mpe_r[++ cnt] = 0;
        poly_mod(mpe_r + mpe_r_begin[tot], a, mpe_mul + mpe_begin[tot], last, len + 1);
        if(r - l + 1 <= shit)
        {
            rep(i, l, r) f[i] = poly_spe(mpe_r + mpe_r_begin[tot], len, b[i]);
            return ;
        }
        if(l == r) return (void)(f[l] = mpe_r[mpe_r_begin[tot]]);
        rep(i, 0, len - 1) a[i] = mpe_r[mpe_r_begin[tot] + i];
        register int mid = (l + r) >> 1;
        poly_mpe_done(f, a, b, len, l, mid, tot << 1);
        rep(i, 0, len - 1) a[i] = mpe_r[mpe_r_begin[tot] + i];
        poly_mpe_done(f, a, b, len, mid + 1, r, tot << 1 | 1);
    }

    int mpe_tmp_a[maxn];



    inline void poly_mpe(int *f, int *a, int *b, int len_a, int len_b)
	{
        register int flag = 0;
        if(!cnt) flag = 1;
        if(flag) poly_dac_mul(b, 1, len_b, 1);
        cnt = 0;
        rep(i, 0, len_a - 1) mpe_tmp_a[i] = a[i];
        poly_mpe_done(f, mpe_tmp_a, b, len_a, 1, len_b, 1), cnt = 0;
        clr(mpe_r);
        if(flag) clr(mpe_mul);
        rep(i, 0, len_a - 1) mpe_tmp_a[i] = 0;
    }

    int mpi_mul[maxn * 10], mpi_tmp[maxn], mpi_l[maxn];

    int mpi_g[maxn], mpi_p[maxn];

    inline void poly_mpi_done(int *f, int *x, int *y, int l, int r, int tot)
	{
        if(l == r) return (void)(f[0] = 1ll * power(mpi_p[l], mod - 2) * y[l] % mod);
        register int mid = (l + r) >> 1, llen = mid - l + 1, rlen = r - mid, len = r - l + 1;
        poly_mpi_done(f, x, y, l, mid, tot << 1);
        rep(i, 0, llen - 1) mpi_mul[i + mpe_begin[tot]] = f[i], f[i] = 0;
        poly_mpi_done(f, x, y, mid + 1, r, tot << 1 | 1);
        rep(i, 0, llen) mpi_tmp[i] = mpe_mul[i + mpe_begin[tot << 1]];
        poly_mul(f, mpi_tmp, len); 
        rep(i, 0, len - 1) mpi_tmp[i] = 0;
        rep(i, 0, rlen) mpi_tmp[i] = mpe_mul[i + mpe_begin[tot << 1 | 1]];
        rep(i, 0, llen - 1) mpi_l[i] = mpi_mul[i + mpe_begin[tot]];
        poly_mul(mpi_l, mpi_tmp, len);
        rep(i, 0, len - 1) mpi_tmp[i] = 0;
        plus(f, mpi_l, len);
        rep(i, 0, len - 1) mpi_l[i] = mpi_mul[i + mpe_begin[tot]] = 0;
    }

    inline void poly_mpi(int *f, int *x, int *y, int len)
	{
        poly_dac_mul(x, 1, len, 1);
        rep(i, 0, len) mpi_g[i] = mpe_mul[i + mpe_begin[1]];
        derivative(mpi_g, len + 1), poly_mpe(mpi_p, mpi_g, x, len, len), poly_mpi_done(f, x, y, 1, len, 1);
        rep(i, 0, len) mpi_g[i] = mpi_p[i] = 0;
    }
}

using polynomial :: inversion;
using polynomial :: poly_pow;
using polynomial :: power;
using polynomial :: init;

int n, m;

int g[maxn], f[maxn];

signed main()
{
    init();
    n = read(), m = read();
    g[0] = 1;
    rep(i, 1, m) g[read() - 1] = mod - 1;
    inversion(f, g, n), poly_pow(f, n, n);
    printf("%lld\n", 1ll * power(n, mod - 2) * f[n - 1] % mod);
    return 0;
}