[多项式ln][多项式exp][多项式求幂][生成函数][DP][FNT] BZOJ 3684: 大朋友和多叉树

$Solution$

把DP写成生成函数的形式。

$$f(x)=x+\sum_{d\in D}f^d(x)$$ 设 $g\big(f(x)\big)=x$，有 $$\begin{eqnarray*} g\big(f(x)\big) &= & f(x)-\sum_{d\in D}f^d(x) \\ g(x) & = &x - \sum_{d\in D}x^d \end{eqnarray*}$$ 根据拉格朗日反演 $$[x^n]h(f(x))={1\over n}[\omega^{n-1}]h'(\omega)({\omega \over g(\omega)})^n$$ 令 $h(x)=x$就得到了 $$[x^n]f(x)={1\over n}[\omega^{n-1}]({\omega\over g(\omega)})^n$$ 现在要求的就是 $({g(\omega)\over \omega})^{-n}$这个东西。 $$f^k(x)=e^{k\ln f(x)}$$ $\ln$直接多项式求逆， $\exp$牛顿迭代。 $\mathcal{O}(n\log n)$。

#include <bits/stdc++.h>
#define show(x) cerr << #x << " = " << x << endl
using namespace std;
typedef long long ll;
typedef pair<int, int> pairs;

const int N = 303030;
const int MOD = 950009857;

inline char get(void) {
    static char buf[100000], *S = buf, *T = buf;
    if (S == T) {
        T = (S = buf) + fread(buf, 1, 100000, stdin);
        if (S == T) return EOF;
    }
    return *S++;
}
template<typename T>
inline void read(T &x) {
    static char c; x = 0; int sgn = 0;
    for (c = get(); c < '0' || c > '9'; c = get()) if (c == '-') sgn = 1;
    for (; c >= '0' && c <= '9'; c = get()) x = x * 10 + c - '0';
    if (sgn) x = -x;
}

inline int pwr(int a, int b) {
    int c = 1;
    while (b) {
        if (b & 1) c = (ll)c * a % MOD;
        b >>= 1; a = (ll)a * a % MOD;
    }
    return c;
}
inline int ivs(int x) {
    return pwr(x, MOD - 2);
}
inline int sum(int a, int b) {
    a += b;
    return a >= MOD ? a - MOD : a;
}
inline int sub(int a, int b) {
    return a < b ? a - b + MOD : a - b;
}
inline void add(int &x, int a) {
    x = sum(x, a);
}

namespace FNT {
    const int MAXN = 303030;
    int ww[MAXN], iw[MAXN];
    int rev[MAXN];
    int num;
    inline void pre(int n) {
        num = n;
        int g = pwr(7, (MOD - 1) / n);
        ww[0] = iw[0] = 1;
        for (int i = 1; i < num; i++)
            iw[n - i] = ww[i] = (ll)ww[i - 1] * g % MOD;
    }
    inline void fnt(int *a, int n, int f) {
        static int x, y, *w;
        w = (f == 1) ? ww : iw;
        for (int i = 0; i < n; i++)
            if (rev[i] > i)
                swap(a[rev[i]], a[i]);
        for (int i = 1; i < n; i <<= 1)
            for (int j = 0; j < n; j += (i << 1))
                for (int k = 0; k < i; k++) {
                    x = a[j + k];
                    y = (ll)a[j + k + i] * w[num / (i << 1) * k] % MOD;
                    a[j + k] = sum(x, y);
                    a[j + k + i] = sub(x, y);
                }
        if (f == -1){
            int in = ivs(n);
            for (int i = 0; i < n; i++)
                a[i] = (ll)a[i] * in % MOD;
        }
    }
}

int inv[N];

inline void pre(int n) {
    inv[1] = 1;
    for (int i = 2; i <= n; i++)
        inv[i] = (ll)(MOD - MOD / i) * inv[MOD % i] % MOD;
}

void getInv(int *a, int *b, int n) {
    using namespace FNT;
    static int tmp[N];
    if (n == 1) return (void)(b[0] = ivs(a[0]));
    getInv(a, b, n >> 1);
    for (int i = 0; i < n; i++) {
        tmp[i] = a[i]; tmp[i + n] = 0;
    }
    int L = 0; while (!(n >> L & 1)) L++;
    for (int i = 0; i < (n << 1); i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L);
    fnt(tmp, n << 1, 1); fnt(b, n << 1, 1);
    for (int i = 0; i < (n << 1); i++)
        tmp[i] = (ll)b[i] * sub(2, (ll)tmp[i] * b[i] % MOD) % MOD;
    fnt(tmp, n << 1, -1);
    for (int i = 0; i < n; i++) {
        b[i] = tmp[i]; b[n + i] = 0;
    }
}
inline void getLn(int *a, int *b, int n) {
    using namespace FNT;
    static int da[N], ia[N], tmp[N];
    for (int i = 0; i < (n << 1); i++)
        tmp[i] = da[i] = ia[i] = 0;
    getInv(a, ia, n);
    for (int i = 1; i < n; i++)
        da[i - 1] = (ll)a[i] * i % MOD;
    int L = 0; while (!(n >> L & 1)) L++;
    for (int i = 0; i < (n << 1); i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L);
    fnt(da, n << 1, 1); fnt(ia, n << 1, 1);
    for (int i = 0; i < (n << 1); i++)
        tmp[i] = (ll)da[i] * ia[i] % MOD;
    fnt(tmp, n << 1, -1);
    b[0] = b[n] = 0;
    for (int i = 1; i < n; i++) {
        b[i] = (ll)tmp[i - 1] * inv[i] % MOD;
        b[n + i] = 0;
    }
}
inline void getExp(int *a, int *b, int n) {
    using namespace FNT;
    static int lb[N];
    if (n == 1) return (void)(b[0] = 1);
    getExp(a, b, n >> 1);
    for (int i = 0; i < (n << 1); i++) lb[i] = 0;
    getLn(b, lb, n);
    for (int i = 0; i < n; i++)
        lb[i] = sub(a[i], lb[i]);
    lb[0] = sum(lb[0], 1);
    int L = 0; while (!(n >> L & 1)) L++;
    for (int i = 0; i < (n << 1); i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << L);
    fnt(b, n << 1, 1); fnt(lb, n << 1, 1);
    for (int i = 0; i < (n << 1); i++)
        b[i] = (ll)b[i] * lb[i] % MOD;
    fnt(b, n << 1, -1);
    for (int i = 0; i < n; i++) b[i + n] = 0;
}
inline void getPow(int *a, int *b, int n, int k) {
    static int la[N];
    k = (k % MOD + MOD) % MOD;
    getLn(a, la, n);
    for (int i = 0; i < n; i++) la[i] = (ll)la[i] * k % MOD;
    getExp(la, b, n);
}

int n, m, x, l;
int g[N], f[N], g_n[N];

int main(void) {
    freopen("1.in", "r", stdin);
    freopen("1.out", "w", stdout);
    FNT::pre(1 << 18);
    pre(1 << 18);
    read(n); read(m);
    for (l = 1; l <= n; l <<= 1);
    for (int i = 0; i < m; i++) {
        read(x); g[x - 1] = MOD - 1;
    }
    g[0] = 1;
    getPow(g, g_n, l, -n);
    cout << (ll)g_n[n - 1] * inv[n] % MOD << endl;
    return 0;
}