生成函数

引入

考虑一类组合对象组成的集合 $A$，对其每个元素定义大小， 记大小为 $n$ 的元素数量为 $A_n$，则 $A$ 的一般生成函数为 $A(x) = \sum_{n\ge 0} A_nx^n$，这里 $A(x$) 是一个形式幂级数，不考虑其收敛域

考虑 $A$ 中元素按顺序组合（笛卡尔积）得到的序列，全部序列所形成的集合 $B$ 满足 $B(x) = \sum_{k\ge0} A(x)^k = \frac{1}{1−A(x)}$

考虑 $A$ 中元素无序组合得到的序列，全部序列所形成的集合 $B$ 满足 $B(x) = \sum_{k\ge0} \frac{A(x)^k}{k!} = \text{ exp }(A(x))$，这里相当于枚举了 $k!$ 种可能的顺序然后去重。

常见生成函数

$$\begin{aligned} \frac{1}{1-x} &= \sum_{k\ge 0} x^k \\ (1+x)^m &= \sum_{k\ge 0} \binom{m}{k}x^k\\ \frac{1}{(1-x)^m} &= \sum_{k\ge 0}\binom{k+m-1}{m-1}x^k(\text{ 考虑隔板法 }) \\ e^x &= \sum_{k\ge 0}\frac{x^k}{k!}\\ e^{-x} &= \sum_{k\ge 0}(-1)^k\frac{x^k}{k!}\\ \ln(1-x) &= -\sum_{k\ge 0} \frac{x^k}{k} \end{aligned}$$

错排的指数生成函数

视一个 $n$ 元排列为带标号的一群置换环的组合，则错排即为没有大小为 1 的置换环。

写出一个环排列的指数生成函数：

$$G(x)=\sum_{n\ge 0} \frac{n!}{n}\times \frac{x^n}{n!}=\sum_{n\ge 0} \frac{x^n}{n} =\ln\frac{1}{1-x}$$

去掉大小为 1 的环，带标号的组合即为

$$F(x)=\text{ exp } (G(x) - x)=\frac{e^{-x}}{1-x}$$

第一类斯特林数

$S(n,k)$ 表示将 $n$ 个数的排列划分成 $k$ 个环的方案数。

考虑一个长度为 $n$ 的环排列的方案数 $c_n=(n-1)!$ ，记其 EGF 为 $C(x)$ , $C(x) = \sum_{n\ge 1} \frac{c_n}{n!}x^n$ 。

对于一个固定的 $k$ ，我们令 $f(n)=S(n,k)$ ，记其 EGF 为 $F(x)$ ， 则有 $F(x)=\frac{C(x)^k}{k!}$ ，同样的我们有 $f(n)=n![x^n]F(x)$。

排列中环个数的期望

沿用上一题的 $C(x)$ , 我们写出 EGF $G(x)$ 表示每个排列的贡献

$$G(x)=\sum_{k\ge 1}\frac{k}{k!}C(x)^k=C(x)\sum_{k\ge 0}\frac{C(x)^k}{k!}=C(x)\text{ exp}(C(x))$$

可以更具体地写出 $C(x)$

$$C(x)=\sum_{k\ge 1} \frac{(k-1)!}{k!} x^k=\sum_{k\ge 1}\frac{1}{k}x^k=-\ln(1 - x)$$

那么可以得到 $G(x) = \frac{-\ln(1-x)}{1-x}$

由于 $[x^n](-\ln(1-x))=\frac{1}{n}$ ，再乘以 $\frac{1}{1-x}$ 实际上是做了一个前缀和。

所以可以得到

$$[x^n]G(x)=\sum_{k\ge 1} \frac{1}{k}$$

有标号简单环森林计数

统计 $n$ 个的无向图，每个点度数都为 $2$ 的方案数。

考虑图中环的计数 , 设 $d_n$ 表示长度为 $n$ 的环的数量，记其 EGF 为 $D(x)$。

$$D(x)=\sum_{k\ge 3} \frac{\frac{(k-1)!}{2}}{k!}x^k=\frac{1}{2}\sum_{k\ge 3} \frac{1}{k}x^k=\frac{-\ln(1-x)-x-\frac{x^2}{2}}{2}$$

记答案的 EGF 为 $G(x)$

$$G(x)=\text{ exp}(D(x))=\frac{e^{-\frac{x}{2}-\frac{x^2}{4}}}{\sqrt{1-x}}$$

有标号简单连通无向图计数

算法1

设 $f(n)$ 表示 $n$ 个点简单无向连通图的个数 ，$g(n)$ 表示 $n$ 个点的图的个数。那么有 $g(n)=2^{\binom{n}{2}}$

枚举第一个点所在的连通块大小，则有：

$$g(n)=\sum_{i=1}^{n}\binom{n-1}{i-1}f(i)g(n-i)$$

移项可得：

$$\frac{g(n)}{(n-1)!}=\sum_{i=1}^{n}\frac{f(i)}{(i-1)!}\times\frac{g(n-i)}{(n-i)!}$$

已知 $g(n)$ 目标求得 $f(n)$

$$\begin{aligned} F(x)&=\sum_{i\ge1} \frac{f(i)}{(i-1)!}x^i\\ G(x)&=\sum_{i\ge0} \frac{g(i)}{i!}x^i\\ H(x)&=\sum_{i\ge1} \frac{g(i)}{(i-1)!}x^i\\ F*G&=H\\ F&=H*G^{-1} \end{aligned}$$

算法2

同样的设法，考虑一张图是由带标号的很多无向连通图组合而成的

设 $f(x)$ , $g(x)$ 分别为 $f$ , $g$ 的指数生成函数

$$\begin{aligned} g(x) &= \text{exp}(f(x))\\ f(x) &= \ln(g(x)) \end{aligned}$$

有标号二分图计数

考虑一个暴力的二分图染色方案数

$$f_n=\sum_{i=0}^{n}\binom{n}{i}2^{(n-i)i}$$

暴力地分成两边，并手动地为二分图定下左右的颜色。

对于这样统计出来的结果，不难发现当一张图中出现 $k$ 个连通块时，总共出现了 $2^k$ 次种连通块染色方案。

对于二分连通图的 EGF $G(x)$ ，设 $\{f_n\}$ 的 EGF 为 $F(x)$。

则有：

$$F(x)=\sum_{k\ge 0}\frac{2^k\times G(x)^k}{k!}=\text{exp }(2G(x))$$

$$\frac{\ln(F(x))}{2}=G(x)$$

又由二分图是二分连通图的组合可得：

$$H(x)=\text{exp}(G(x))=\sqrt{F(x)}$$

另一个视角

对于 $F(x)$ 我们取 $\ln$ 可以求得 带颜色的二分图连通图数量的 EGF ，我们记为 $B(x)$。

发现二分图本身是没有颜色的，因此一个二分图会被记成两次 所以有 $G(x)=\frac{B(x)}{2}$

同样的有

$$H(x)=\text{exp}(G(x))=\sqrt{F(x)}$$

分拆数

总共有 $n$ 种物品，每种物品有 $i$ 的体积，求填满大小为 $n$ 的背包的方案数。

解法1

写出生成函数

$$\begin{aligned} P(x)&=\prod_{i=1}^n \sum_{j\ge 0}x^{j\times i} \\ &=\prod_{i=1}^n \frac{1}{1-x^i} \end{aligned}$$

取对数可得：

$$\begin{aligned} \ln P(x)&=-\sum_{i=1}^n \ln(1 -x^i)\\ &=\sum_{i=1}^n\sum_{j\ge 0}\frac{x^{j\times i}}{j}\\ &=\sum_{i=1}^n x^i \sum_{j | i} \frac{1}{j} \end{aligned}$$

根据调和级数，可以通过 $O(n\log n)$ 的复杂度求出 $\ln P(x)$ 然后 $\text{exp}$ 回去即可得到 $P(x)$.

一道组合恒等式

给定 $n$ 求：

$$\sum_{k\ge 0}^{n} \binom{k}{n-k}$$

记答案为 $f(n)$，写出 OGF $F(x)$ :

$$\begin{aligned} F(x) &= \sum_{n\ge 0} \sum_{k\ge 0}^{n} \binom{k}{n - k}x^n=\sum_{k\ge 0}\sum_{n\ge k} \binom{k}{n - k}x^n\\ &=\sum_{k\ge 0}x^k\sum_{n-k\ge 0}\binom{k}{n - k}x^{n-k}\\ &=\sum_{k\ge 0}x^k(1+x)^k\\ &=\frac{1}{1-x-x^2} \end{aligned}$$

发现最后的形式是斐波那契数列的 OGF

[TJOI2015]概率论

设 $g(i)$ 表示 $i$ 个结点的二叉树总数，$f(i)$ 表示 $i$ 个结点的二叉树的叶子结点数的和。

目标求得 $f(n)$ 时考虑每个 $f(i)$ 的贡献，枚举左子树大小，根据对称性，最后乘2即可得到所有的贡献。

$$f(n)=2\sum_{i=1}^{n-1}f(i)g(n-1-i)$$

写成生成函数形式可得：

$$F(x)=\sum_{n\ge 0}f(n)x^n=0+x+\sum_{n\ge0}2x\sum_{i=1}^{n-1}f(i)x^i\times g(n-1-i)x^{n-1-i}$$

$$F(x)=x+2xF(x)G(x)$$

已知 $G(x)$ 为卡特兰数的 OGF，由其封闭形式 $G(x)=\frac{1-\sqrt{1-4x}}{2x}$得

$$F(x)=\frac{x}{\sqrt{1-4x}}$$

暴力展开得：

$$\begin{aligned} F(x)&=x(1+(-4x))^{-\frac{1}{2}}=x\sum_{k\ge0}\binom{-\frac{1}{2}}{k}(-4x)^k\\ &=x\sum_{k\ge0}\frac{\prod_{i=0}^{k-1}(-\frac{1}{2}-i)}{k!}4^k(-x)^k\\ &=x\sum_{k\ge0}\frac{(2k-1)!!}{k!}2^kx^k\\ &=x\sum_{k\ge 0}\frac{(2k)!}{k!k!}x^k\\ &=\sum_{k\ge1}\binom{2k-2}{k-1}x^k \end{aligned}$$

最后求的是期望所以是：

$$\frac{[x^n]F(x)}{[x^n]G(x)}=\frac{\binom{2n-2}{n-1}}{\frac{\binom{2n}{n}}{n+1}}=\frac{n^2+n}{4n-2}$$

所以有 $f(n)=f_{n+1}$ 最后答案为斐波那契数列的第 $n+1$ 项。

CF891 E

你有 $n$ 个数 $a_1,a_2,\dots,a_n$ 要进行 $k$ 次操作，每次随机选择一个数 $x \in [1,n]$，把 $a_x$ 减一，并将答案增加除 $a_x$ 外所有数的乘积。

求最终答案的期望，答案对 $10^9+7$ 取模。

列出一次操作的贡献：

$$\prod_{i=1,i\neq x}^na_i=\prod_{i=1}^n a_i-(a_x-1)\times\prod_{i=1,i\neq x}^na_i$$

再往后写几项相加发现中间项都被消掉了，设每个数 $a_i$ 在 $k$ 轮操作后被减少的数值为 $b_i$。

那么对于一种合法局面的贡献值为：

$$\prod_{i=1}^n a_i-\prod_{i=1}^n (a_i - b_i)$$

考虑出现的概率，首先是 $\sum b_i = k$ 的方案数是可重集排列数 $\binom{k}{b_1,b_2,\dots,b_n}$，总方案数是每次任意挑一个数减 $n^k$

列出期望的式子：

$$\Epsilon = \prod_{i=1}^n a_i - \frac{k!}{n^k}\prod_{i=1}^n\frac{(a_i-b_i)}{b_i!}$$

我们关心的是 $b_i$ 的具体值，需要满足 $\sum b_i=k$，可以枚举 $b_i$ 列出 EGF。

$$\begin{aligned} \Epsilon &= \prod_{i=1}^n a_i - \frac{k!}{n^k}\prod_{i=1}^n\sum_{j\ge 0}\frac{(a_i-j)}{j!}\times x^j\\ &= \prod_{i=1}^n a_i - \frac{k!}{n^k}\prod_{i=1}^n(a_i-x)\times e^x\\ &= \prod_{i=1}^n a_i - \frac{k!}{n^k}\times e^{nx}\prod_{i=1}^n(a_i-x) \end{aligned}$$

实际上我们只关心 $[x^k]\Epsilon$，上面式子复杂度的瓶颈在于后半部分，可以分治FFT $O(n\log n)$ 求出 。

AGC005F Many Easy Problems

题意

给出一个 $n$ 个点的树，对每一个 $k$ 求出 $f(k)$。

$f(k)$ 为所有大小为 $k$ 的点集 $S$ ，能覆盖 $S$ 的最小连通块大小的和。

$$n\le 200000$$

题解

不能够直接按联通块枚举。

考虑求问题的补，即求每一个点不在多少种选取 $S$ 的方案中做贡献。

假设当前考虑的一个点 $u$ 为根，那么当且仅当选取的 $k$ 个点都在 $u$ 同一个儿子的子树中，那么 $u$ 不做贡献的方案数为：

$$\sum_{v\in \text{son}_u} \binom{\text{sz}_v}{k}$$

可以得到

$$f(k)=n\binom{n}{k}-\sum_{u=1}^{n}\sum_{v\in \text{son}_u} (\binom{\text{sz}_v}{k}+\binom{n-\text{sz}_v}{k})$$

实际上我们现在不需要关心树的形态，发现每个 $\text{sz}_v$ 仅会被他的父亲贡献。那么现在我们就只需要考虑每种 $\text{sz}$ 的出现次数了。

$$f(k)=n\binom{n}{k}-\sum_{i=1}^{n}\text{cnt}_i\binom{i}{k}$$

现在目标则是对于每个 $k$ 求出 $\sum_{i=1}^{n}\text{cnt}_i\binom{i}{k}$。写出其 OGF:

$$\begin{aligned} F(x)&=\sum_{k=1}^{n}\sum_{i=1}^{n}\text{cnt}_i\binom{i}{k}x^k\\ &=\sum_{k=1}^{n}\frac{x^k}{k!}\sum_{i=k}^{n}(i!\text{ cnt}_i) \times \frac{1}{(i-k)!} \end{aligned}$$

发现后面的是一个差卷积的形式，考虑构造一个 OGF 的转置

记 $A(x)=\sum_{i=0}^n (i!\text{ cnt}_i) x^i$ , $B(x)=\sum_{i=0}^{n}\frac{1}{(n-i)!}x^i$

令 $H(x)= A(x)\times B(x)$，那么 $[x^{k+n}]H(x)$ 即为我们所要求的 $k![x^k]F(x)$

代码

#include <bits/stdc++.h>

using namespace std;

using int64 = long long;

template<int MOD>
struct ModInt {
    int x;
 
    ModInt() : x(0) {}
    ModInt(int64 y) : x(y >= 0 ? y % MOD : (MOD - (-y) % MOD) % MOD) {}

    inline int inc(const int &x) {
        return x >= MOD ? x - MOD : x;
    }
    inline int dec(const int &x) {
        return x < 0 ? x + MOD : x;
    }
 
    ModInt &operator+= (const ModInt &p) {
        x = inc(x + p.x);
        return *this;
    } 
    ModInt &operator-= (const ModInt &p) {
        x = dec(x - p.x);
        return *this;
    }
 
    ModInt &operator*= (const ModInt &p) {
        x = (int)(1ll * x * p.x % MOD);
        return *this;
    }
    ModInt &operator/= (const ModInt &p) {
        *this *= p.inverse();
        return *this;
    }
 
    ModInt operator-() const { return ModInt(-x); } 

    friend ModInt operator + (const ModInt& lhs, const ModInt& rhs) {
        return ModInt(lhs) += rhs;
    }
    friend ModInt operator - (const ModInt& lhs, const ModInt& rhs) {
        return ModInt(lhs) -= rhs;
    }
    friend ModInt operator * (const ModInt& lhs, const ModInt& rhs) {
        return ModInt(lhs) *= rhs;
    }
    friend ModInt operator / (const ModInt& lhs, const ModInt& rhs) {
        return ModInt(lhs) /= rhs;
    }
 
    bool operator == (const ModInt &p) const { return x == p.x; } 
    bool operator != (const ModInt &p) const { return x != p.x; }
 
    ModInt inverse() const {
        int a = x, b = MOD, u = 1, v = 0, t;
        while(b > 0) {
            t = a / b;
            swap(a -= t * b, b);
            swap(u -= t * v, v);
        }
        return ModInt(u);
    }

    ModInt pow(int64 n) const {
        ModInt ret(1), mul(x);
        while(n > 0) {
            if(n & 1) ret *= mul;
            mul *= mul;
            n >>= 1;
        }
        return ret;
    }
 
    friend ostream &operator<<(ostream &os, const ModInt &p) {
        return os << p.x;
    }
 
    friend istream &operator>>(istream &is, ModInt &a) {
        int64 t;
        is >> t;
        a = ModInt<MOD>(t);
        return (is);
    }
    static int get_mod() { return MOD; }
};

const int MOD = 924844033;

using MODint = ModInt<MOD>;
MODint pow(int64 n, int64 x) { return MODint(n).pow(x); }
MODint pow(MODint n, int64 x) { return n.pow(x); }

vector<int> rev;
vector<MODint> roots{0 , 1};

void dft(vector<MODint> &a) {
    int n = a.size();

    if (rev.size() != n ) {
        int k = __builtin_ctz(n) - 1;
        rev.resize(n);
        for (int i = 0 ; i < n ; i++) {
            rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
        }
    }

    for (int i = 0 ; i < n ; i++) {
        if (rev[i] < i) {
            swap(a[i] , a[rev[i]]);
        }
    }

    if (roots.size() < n) {
        int k = __builtin_ctz(roots.size());
        roots.resize(n);
        while ((1 << k) < n) {
            MODint e = pow(MODint(5) , (MOD - 1) >> (k + 1));
            for (int i = 1 << (k - 1) ; i < (1 << k) ; i++) {
                roots[2 * i] = roots[i];
                roots[2 * i + 1] = roots[i] * e;
            }
            ++k;
        }
    }

    for (int k = 1 ; k < n ; k <<= 1) {
        for (int i = 0 ; i < n ; i += 2 * k) {
            for (int j = 0 ; j < k ; j++) {
                MODint u = a[i + j];
                MODint v = a[i + j + k] * roots[k + j];
                a[i + j] = u + v;
                a[i + j + k] = u - v;
            }
        }
    }
}

void idft(vector<MODint> &a) {
    int n = a.size();
    reverse(a.begin() + 1 , a.end());
    dft(a);
    MODint inv = (1 - MOD) / MODint(n);
    for (int i = 0 ; i < n ; i++) {
        a[i] *= inv;
    }
}

struct poly {
    vector<MODint> a;
    poly() {} 
    poly(const vector<MODint> &a) : a(a) {}

    int size() const {
        return a.size();
    }
    void resize(int n) {
        a.resize(n);
    }
    MODint operator[](int idx) const {
        if (idx < 0 || idx >= size()) {
            return 0;
        }
        return a[idx];
    }
    MODint &operator[](int idx) {
        return a[idx];
    }

    poly mulxk(int k) const {
        auto b = a;
        b.insert(b.begin(), k, 0);
        return poly(b);
    }
    poly modxk(int k) const {
        k = min(k, size());
        return poly(vector<MODint>(a.begin(), a.begin() + k));
    }
    poly divxk(int k) const {
        if (size() <= k) {
            return poly();
        }
        return poly(vector<MODint>(a.begin() + k, a.end()));
    }

    friend poly operator+(const poly &a, const poly &b) {
        vector<MODint> res(max(a.size(), b.size()));
        for (int i = 0; i < int(res.size()); i++) {
            res[i] = a[i] + b[i];
        }
        return poly(res);
    }
    friend poly operator-(const poly &a, const poly &b) {
        vector<MODint> res(max(a.size(), b.size()));
        for (int i = 0; i < int(res.size()); i++) {
            res[i] = a[i] - b[i];
        }
        return poly(res);
    }
    friend poly operator*(poly a, poly b) {
        if (a.size() == 0 || b.size() == 0) {
            return poly();
        }
        int sz = 1, tot = a.size() + b.size() - 1;
        while (sz < tot)
            sz *= 2;
        a.a.resize(sz);
        b.a.resize(sz);
        dft(a.a);
        dft(b.a);
        for (int i = 0; i < sz; ++i) {
            a.a[i] = a[i] * b[i];
        }
        idft(a.a);
        a.resize(tot);
        return a;
    }
    friend poly operator*(MODint a, poly b) {
        for (int i = 0; i < int(b.size()); i++) {
            b[i] *= a;
        }
        return b;
    }
    friend poly operator*(poly a, MODint b) {
        for (int i = 0; i < int(a.size()); i++) {
            a[i] *= b;
        }
        return a;
    }
    poly &operator+=(poly b) {
        return (*this) = (*this) + b;
    }
    poly &operator-=(poly b) {
        return (*this) = (*this) - b;
    }
    poly &operator*=(poly b) {
        return (*this) = (*this) * b;
    }
    poly deriv() const {
        if (a.empty()) {
            return poly();
        }
        vector<MODint> res(size() - 1);
        for (int i = 0; i < size() - 1; ++i) {
            res[i] = (i + 1) * a[i + 1];
        }
        return poly(res);
    }
    poly integr() const {
        vector<MODint> res(size() + 1);
        for (int i = 0; i < size(); ++i) {
            res[i + 1] = a[i] / (i + 1);
        }
        return poly(res);
    }
    poly inv(int m) const {
        poly x({a[0].inverse()});
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (poly({2}) - modxk(k) * x)).modxk(k);
        }
        return x.modxk(m);
    }
    poly log(int m) const {
        return (deriv() * inv(m)).integr().modxk(m);
    }
    poly exp(int m) const {
        poly x({1});
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (poly({1}) - x.log(k) + modxk(k))).modxk(k);
        }
        return x.modxk(m);
    }
    poly sqrt(int m) const {
        poly x({1});
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x + (modxk(k) * x.inv(k)).modxk(k)) * ((MOD + 1) / 2);
        }
        return x.modxk(m);
    }
    poly mulT(poly b) const {
        if (b.size() == 0) {
            return poly();
        }
        int n = b.size();
        reverse(b.a.begin(), b.a.end());
        return ((*this) * b).divxk(n - 1);
    }
    vector<MODint> eval(vector<MODint> x) const {
        if (size() == 0) {
            return vector<MODint>(x.size(), 0);
        }
        const int n = max(int(x.size()), size());
        vector<poly> q(4 * n);
        vector<MODint> ans(x.size());
        x.resize(n);
        function<void(int, int, int)> build = [&](int p, int l, int r) {
            if (r - l == 1) {
                q[p] = poly({1, -x[l]});
            } else {
                int m = (l + r) / 2;
                build(2 * p, l, m);
                build(2 * p + 1, m, r);
                q[p] = q[2 * p] * q[2 * p + 1];
            }
        };
        build(1, 0, n);
        function<void(int, int, int, const poly &)> work = [&](int p, int l, int r, const poly &num) {
            if (r - l == 1) {
                if (l < int(ans.size())) {
                    ans[l] = num[0];
                }
            } else {
                int m = (l + r) / 2;
                work(2 * p, l, m, num.mulT(q[2 * p + 1]).modxk(m - l));
                work(2 * p + 1, m, r, num.mulT(q[2 * p]).modxk(r - m));
            }
        };
        work(1, 0, n, mulT(q[1].inv(n)));
        return ans;
    }
};

template<int SIZE , int MOD>
struct Combine {
    ModInt<MOD> fac[SIZE + 5] , inv[SIZE + 5];
    Combine() {
        fac[0] = 1;
        for (int i = 1 ; i <= SIZE ; i++) {
            fac[i] = fac[i - 1] * i;
        }
        inv[SIZE] = fac[SIZE].inverse();
        for (int i = SIZE ; i >= 1 ; i--) {
            inv[i - 1] = inv[i] * i;
        }
    }
    inline ModInt<MOD> operator()(int n , int m) {
        if (n < m) return 0;
        return fac[n] * inv[m] * inv[n - m];
    }
};

Combine<200000 , MOD> C;

inline MODint fac(int n) { return C.fac[n]; }
inline MODint inv(int n) { return C.inv[n]; }

int main() {
    cin.tie(nullptr)->sync_with_stdio(0);

    int n;
    cin >> n;

    vector<vector<int>> g(n + 1 , vector<int>());

    for (int i = 1 ; i < n ; i++) {
        int u , v;
        cin >> u >> v;
        g[u].push_back(v);
        g[v].push_back(u);
    } 

    vector<int> sz(n + 1) , cnt(n + 1);

    function<void(int , int)> dfs = [&](int u , int fa){
        sz[u] = 1;
        for (int v : g[u]) {
            if (v == fa) continue;
            dfs(v , u);
            ++cnt[sz[v]];
            ++cnt[n - sz[v]];
            sz[u] += sz[v];
        }
    };

    dfs(1 , 0);

    poly a , b;
    a.resize(n + 1);
    b.resize(n + 1);

    for (int i = 0 ; i <= n ; i++) {
        a[i] = cnt[i] * fac(i);
    }
    for (int i = 0 ; i <= n ; i++) {
        b[i] = inv(n - i);
    }
    poly c = a * b;

    for (int i = 1 ; i <= n ; i++) {
        cout << n * C(n , i) - c[i + n] * inv(i) << '\n';
    }
    return 0;
}

ICPC Shenyang 2018 M Renaissance Past in Nancy

题意

有 $n$ 个物品，每个物品有数量 $a_i$ 和 价格 $b_i$，每次询问一个区间 $[l,r]$ ，和当前的钱的数量 $c$ ，求区间的一个背包方案数，强制在线。

$$n,m\le 10000\;a_i,b_i\le 1000 \; c\le1000$$

题解

不难写出一个物品的 OGF 为：

$$\frac{1-x^{(a_i+1)\times b_i}}{1-x^b}$$

对于最终的答案，可以列出下式：

$$\text{ans} =\prod_{i=l}^r \frac{1-x^{(a_i+1)\times b_i}}{1-x^{b_i}}$$

可以预处理出 $F_i(x)$ 和 $F^{-1}_i(x)$:

$$\begin{aligned} F_i(x)&=\prod_{j=1}^i \frac{1-x^{(a_i+1)\times b_i}}{1-x^{b_i}}\\ F^{-1}_i(x)&=\prod_{j=1}^i \frac{1-x^{b_i}}{1-x^{(a_i+1)\times b_i}}\\ \end{aligned}$$

看起来只要暴力乘出来就好了，考虑怎么实现这个过程。

发现分子的乘法转移就是一个移位相减，分母似乎不大好处理。

回归 OGF $\frac{1}{1-x^{b_i}}$ 的本质：

$$\frac{1}{1-x^{b_i}}=\sum_{j\ge 0}x^{j\times b_i}$$

考虑背包问题，一个无限的物品，又回归到完全背包上。因此对于分母的转移我们只需类似完全背包转移即可。

注意我们所求得是：

$$\sum_{i=0}^c[x^i]F_r(x)\times F_{l-1}^{-1}(x)$$

对 $F_i(x)$ 做一次前缀和即可。

PS：做一个前缀和相当于乘上一个 $\frac{1}{1-x}$

代码

int main() {
    cin.tie(nullptr);

    int T;
    cin >> T;

    for (int CASE = 1 ; CASE <= T ; CASE++) {
        cout << "Case #" << CASE << ":\n";

        int n , m;
        cin >> n >> m;

        vector<vector<MODint>> f(n + 1 , vector<MODint>(1001));
        vector<vector<MODint>> g(n + 1 , vector<MODint>(1001));

        f[0][0] = g[0][0] = 1;

        for (int i = 1 ; i <= n ; i++) {
            int a , b;
            cin >> a >> b;

            f[i] = f[i - 1];
            g[i] = g[i - 1];

            for (int j = b ; j <= 1000 ; j++) {
                f[i][j] += f[i][j - b];
            }
            for (int j = b * (a + 1) ; j <= 1000 ; j++) {
                g[i][j] += g[i][j - b * (a + 1)];
            }

            for (int j = 1000 ; j >= b * (a + 1) ; j--) {
                f[i][j] -= f[i][j - b * (a + 1)];
            }
            for (int j = 1000 ; j >= b ; j--) {
                g[i][j] -= g[i][j - b];
            }
        }

        // 乘上 1/1-x 相当于做一个重量为1的完全背包
        for (int i = 1 ; i <= n ; i++) {
            for (int j = 1 ; j <= 1000 ; j++) {
                f[i][j] += f[i][j - 1];
            }
        }

        MODint lastans = 0;
        while (m--) {
            int l , r , c , d;
            cin >> l >> r >> c;

            d = lastans.x;
            l = (l + d) % n + 1;
            r = (r + d) % n + 1;
            if (l > r) swap(l , r);

            lastans = 0;
            for (int i = 0 ; i <= c ; i++) {
                lastans += f[r][i] * g[l - 1][c - i];
            }
            cout << lastans << '\n';

        }
    }
    return 0;
}

[SDOI2017] 遗忘的集合

题解

一个集合中出现的数做完全背包得到了 $f(i)$ 的值。

不妨设 $a_i \in [0,1]$ 表示 $i$ 在集合中是否出现。

所以可以列出 $f(i)$ 的生成函数：

$$F(x) = \prod_{i=1}^{n}\frac{1}{(1-x^i)^{a_i}}$$

通过上面提到的分拆数的类似做法，可以得到：

$$\ln F(x) = \sum_{i=1} x^{i}\sum_{j|i} \frac{j\times a_j}{i}$$

$$\ln F(x) = \sum_{i = 1}\frac{x^i}{i}\sum_{j|i}j\times a_j$$

后面是一个经典的反演形式，因为我们已知 $\ln F(x)$ ，求 $g(x)= x\times a_x$，因为 $f = g * I$ 所以有 $g = f * \mu$。

注意多项式求逆的时候，任意模数的取模处理。

每一次做完乘法之后都要重新 CRT 合并一下 ，尤其要注意取负数是要在给定模数的值域下取，而非 NTT 里的模数。除了乘法部分尽可能多地去用原来的模数去取模。

实际上我更倾向于拆系数 FFT

代码

#include <bits/stdc++.h>

using namespace std;

using int64 = long long;

int mpow(int x , int y , int MOD) {
  int res = 1;
  while(y) {
    if(y & 1) res = 1ll * res * x % MOD;
    x = 1ll * x * x % MOD;
    y >>= 1;
  }
  return res;
}


int inverse(long long x , const int &MOD) {
  x %= MOD;
    return mpow(x , MOD - 2 , MOD);
}

const int MOD1 = 998244353 , MOD2 = 1004535809 , MOD3 = 469762049;
const int64 MOD4 = 1ll * MOD1 * MOD2;

const int INV1 = inverse(MOD1 , MOD2);
const int INV2 = inverse(MOD4 % MOD3 , MOD3);

int P;

struct MInt {
    int v1 , v2 , v3;
    
    MInt(int64 y = 0) : v1(y % MOD1) , v2(y % MOD2) , v3(y % MOD3) {}
    MInt(int64 a , int64 b , int64 c) : v1(a % MOD1) , v2(b % MOD2) , v3(c % MOD3) {}

    int inc(const int &x , const int &MOD) const {
        return x >= MOD ? x - MOD : x;
    }
    int dec(const int &x , const int &MOD) const {
        return x < 0 ? x + MOD : x;
    }

    void add(int &x , const int &y , const int &MOD) {
        x = inc(x + y , MOD);
    }
    void sub(int &x , const int &y , const int &MOD) {
        x = dec(x - y , MOD);
    }

    MInt &operator+= (const MInt &p) {
        add(v1 , p.v1 , MOD1);
        add(v2 , p.v2 , MOD2);
        add(v3 , p.v3 , MOD3);
        return *this;
    } 
    MInt &operator-= (const MInt &p) {
        sub(v1 , p.v1 , MOD1);
        sub(v2 , p.v2 , MOD2);
        sub(v3 , p.v3 , MOD3);
        return *this;
    }
 
    MInt &operator*= (const MInt &p) {
        v1 = 1ll * v1 * p.v1 % MOD1;
        v2 = 1ll * v2 * p.v2 % MOD2;
        v3 = 1ll * v3 * p.v3 % MOD3;
        return *this;
    }
    MInt &operator/= (const MInt &p) {
        v1 = 1ll * v1 * ::inverse(p.v1 , MOD1) % MOD1;
        v2 = 1ll * v2 * ::inverse(p.v2 , MOD2) % MOD2;
        v3 = 1ll * v3 * ::inverse(p.v3 , MOD3) % MOD3;
        return *this;
    }

    friend MInt operator + (const MInt& lhs, const MInt& rhs) {
        return MInt(lhs) += rhs;
    }
    friend MInt operator - (const MInt& lhs, const MInt& rhs) {
        return MInt(lhs) -= rhs;
    }
    friend MInt operator * (const MInt& lhs, const MInt& rhs) {
        return MInt(lhs) *= rhs;
    }
    friend MInt operator / (const MInt& lhs, const MInt& rhs) {
        return MInt(lhs) /= rhs;
    }

    MInt unit_pow(int k) {
        return MInt(mpow(3 , (MOD1 - 1) >> (k + 1) , MOD1) ,
                    mpow(3 , (MOD2 - 1) >> (k + 1) , MOD2) ,
                    mpow(3 , (MOD3 - 1) >> (k + 1) , MOD3));
    }

    MInt inverse() const {
        return MInt(::inverse(v1 , MOD1) , ::inverse(v2 , MOD2) , ::inverse(v3 , MOD3));
    }

    MInt idft(int n) {
        return MInt(n).inverse();
    }

    MInt operator-() const {
        return MInt(dec(-v1 , MOD1) , dec(-v2 , MOD2) , dec(-v3 , MOD3));
    }

    int get_val(int P) const {
        int64 x = 1ll * (v2 - v1 + MOD2) % MOD2 * INV1 % MOD2 * MOD1 + v1;
        int64 y = 1ll * (v3 - x % MOD3 + MOD3) % MOD3 * INV2 % MOD3 * (MOD4 % P) % P + x;
        return y % P;
    }
};
vector<int> rev;
vector<MInt> roots{0 , 1};

void dft(vector<MInt> &a) {
    int n = a.size();

    if (rev.size() != n ) {
        int k = __builtin_ctz(n) - 1;
        rev.resize(n);
        for (int i = 0 ; i < n ; i++) {
            rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
        }
    }

    for (int i = 0 ; i < n ; i++) {
        if (rev[i] < i) {
            swap(a[i] , a[rev[i]]);
        }
    }

    if (roots.size() < n) {
        int k = __builtin_ctz(roots.size());
        roots.resize(n);
        while ((1 << k) < n) {
            MInt e = MInt().unit_pow(k);
            for (int i = 1 << (k - 1) ; i < (1 << k) ; i++) {
                roots[2 * i] = roots[i];
                roots[2 * i + 1] = roots[i] * e;
            }
            ++k;
        }
    }

    for (int k = 1 ; k < n ; k <<= 1) {
        for (int i = 0 ; i < n ; i += 2 * k) {
            for (int j = 0 ; j < k ; j++) {
                MInt u = a[i + j];
                MInt v = a[i + j + k] * roots[k + j];
                a[i + j] = u + v;
                a[i + j + k] = u - v;
            }
        }
    }
}

void idft(vector<MInt> &a) {
    int n = a.size();
    reverse(a.begin() + 1 , a.end());
    dft(a);
    MInt inv = (1) / MInt(n);
    for (int i = 0 ; i < n ; i++) {
        a[i] *= inv;
    }
}

struct poly {
    vector<MInt> a;
    poly() {} 
    poly(const vector<MInt> &a) : a(a) {}

    int size() const {
        return a.size();
    }
    void resize(int n) {
        a.resize(n);
    }
    MInt operator[](int idx) const {
        if (idx < 0 || idx >= size()) {
            return 0;
        }
        return a[idx];
    }
    MInt &operator[](int idx) {
        return a[idx];
    }

    poly mulxk(int k) const {
        auto b = a;
        b.insert(b.begin(), k, 0);
        return poly(b);
    }
    poly modxk(int k) const {
        k = min(k, size());
        return poly(vector<MInt>(a.begin(), a.begin() + k));
    }
    poly divxk(int k) const {
        if (size() <= k) {
            return poly();
        }
        return poly(vector<MInt>(a.begin() + k, a.end()));
    }

    friend poly operator+(const poly &a, const poly &b) {
        vector<MInt> res(max(a.size(), b.size()));
        for (int i = 0; i < int(res.size()); i++) {
            res[i] = a[i] + b[i];
        }
        return poly(res);
    }
    friend poly operator-(const poly &a, const poly &b) {
        vector<MInt> res(max(a.size(), b.size()));
        for (int i = 0; i < int(res.size()); i++) {
            res[i] = a[i] - b[i];
        }
        return poly(res);
    }
    friend poly operator*(poly a, poly b) {
        if (a.size() == 0 || b.size() == 0) {
            return poly();
        }
        int sz = 1, tot = a.size() + b.size() - 1;
        while (sz < tot)
            sz *= 2;
        a.a.resize(sz);
        b.a.resize(sz);
        dft(a.a);
        dft(b.a);
        for (int i = 0; i < sz; ++i) {
            a.a[i] = a[i] * b[i];
        }
        idft(a.a);
        a.resize(tot);
        return a;
    }

    poly &operator+=(poly b) {
        return (*this) = (*this) + b;
    }
    poly &operator-=(poly b) {
        return (*this) = (*this) - b;
    }
    poly &operator*=(poly b) {
        return (*this) = (*this) * b;
    }
    poly deriv() const {
        if (a.empty()) {
            return poly();
        }
        vector<MInt> res(size() - 1);
        for (int i = 0; i < size() - 1; ++i) {
            res[i] = MInt(1ll * (i + 1) * a[i + 1].get_val(P) % P);
    }
        return poly(res);
    }
    poly integr() const {
        vector<MInt> res(size() + 1);
        for (int i = 0; i < size(); ++i) {
            res[i + 1] = MInt(1ll * a[i].get_val(P) * ::inverse(i + 1 , P) % P);
        }
        return poly(res);
    }
    poly inv(int m) const {
        poly x({::inverse(a[0].v1 , P)});
        int k = 1;
        while (k < m) {
            k *= 2;
            poly y = modxk(k) * x;
            for (int i = 1 ; i < y.size() ; i++) {
                y[i] = MInt((P - y[i].get_val(P)) % P);
            }
            y[0] = MInt(1);
            x = (y * x).modxk(k);
            for (int i = 0 ; i < x.size() ; i++) {
                x[i] = MInt(x[i].get_val(P));
            }
        }
        return x.modxk(m);
    }
    poly log(int m) const {
        poly p = (deriv() * inv(m)).integr();
    return p.modxk(m);
    }
    poly mulT(poly b) const {
        if (b.size() == 0) {
            return poly();
        }
        int n = b.size();
        reverse(b.a.begin(), b.a.end());
        return ((*this) * b).divxk(n - 1);
    }
};

int main() {
    cin.tie(nullptr)->sync_with_stdio(0);

    int n , l = 1;
    cin >> n >> P;

  while (l <= n) l <<= 1;

    vector<MInt> a(n + 1);

    a[0] = MInt(1);
    for (int i = 1 ; i <= n ; i++) {
        int x;
        cin >> x;
        a[i] = MInt(x);
    }

    poly pa(a);
    pa = pa.log(l);

    vector<int> vis(n + 1) , mu(n + 1) , pri;

    mu[1] = 1;
    for (int i = 2 ; i <= n ; i++) {
        if (!vis[i]) {
            mu[i] = -1;
            pri.push_back(i);
        }
        for (int p : pri) {
            if (1ll * p * i > n) break;
            vis[i * p] = 1;
            if (i % p == 0) {
                mu[i * p] = 0;
                break;
            }
            mu[i * p] = -mu[i];
        }
    }

    vector<int> f(n + 1) , c(n + 1);
    for (int i = 1 ; i <= n ; i++) {
        f[i] = 1ll * pa[i].get_val(P) * i % P;
    }

    for (int i = 1 ; i <= n ; i++) {
        for (int j = i ; j <= n ; j += i) {
            (c[j] += 1ll * f[i] * mu[j / i]) %= P;
            (c[j] += P) %= P;
        }
    }

    int cnt = count_if(c.begin() ,c.end(), [](int x) { return x != 0; });
    cout << cnt << '\n';
    for (int i = 1 ; i <= n ; i++) {
        if (c[i] != 0) {
      cout << i << ' ';
    }
    }

    return 0;
}

牛客-wyf的超级多项式

给出一个 $F_i$ 的计算方式：

$$F_i=\sum_{j=1}^{k}a_j\times v_j^i$$

现在给出 $v_1,v_2,\cdots,v_k$ 和 $F_1,F_2,\cdots,F_k$，求 $F_n$

$$n,k\le 10^5, n-k\le 1000$$

题解

写出 $F_i$ 的生成函数 $F(x)$:

$$F(x)=\sum_{i=1}^{k} \frac{a_i}{1-v_ix}$$

通分可得 $F(x) = \frac{A(x)}{B(x)}$ ，我们不知道 $A(x)$ 实际上我们也不关心 $A(x)$ 的具体形式。可以首先分治 NTT 求出 $B(x)$

令 $B(x) = 1-\sum_{i=1}^{k} c_i x^i$ ，我们记 $C(x)=\sum_{i=1}^k c_i x^i$

于是我们有

$$\begin{aligned} F(x)&=\frac{A(x)}{1-C(x)}\\ F(x)&=A(x)+C(x)F(x) \end{aligned}$$

不难发现 $A(x)$ 的次数小于 $k$ ，那么对于 $\forall n \gt k$ , 都可以通过 $F$ 和 $C$ 递推得到当前的 $F$。

China-Final 2016 I Cherry Pick

题意

总共有 $n$ 颗樱桃树，在每个樱桃树成功摘取的概率为 $p$，每个樱桃树之间独立。

现在手上有 $M$ 种硬币，每种硬币的面值为 $c_i$ ，每种硬币都有无限多个。若摘完 $n$ 颗樱桃树之后，总共取得了 $k$ 个樱桃，那么所要支付的值 $f(k)$ 为最小的大于等于 $k$ 的数 $x$ 减去 $k$。即最少的多花的钱。

求 最少多花的钱的期望。

$$n\le 10^9,M \le 100,c_i\le 10000$$

$$\forall i \neq j,c_i\times c_j\le 10000$$

任意的 $c_i$ 不重复。数据里有不少 $M=1$ 的情况。

题解

这题实际上是在求：

$$\sum_{k\ge 0} \binom{n}{k} p^k(1-p)^{n-k}\times f(k)$$

考虑先求 $f(k)$ ，经典的同余最短路，可以看看 OIWIKI。取最小的 $c_1$ 作为模数，其余的 $c_i$ 作为一条转移边，求得 $d_i$ 表示最小的 $x$ 满足 $x\equiv i\text{ mod } c_1)$。然后可以找到循环节为 $c_1$

考虑前面那个东西怎么求，可以列出个生成函数 $(px+(1-p))^n$ ，发现答案在一定范围内只与 模 $c_i$ 意义下的值有关，所以只要用倍增+FFT求出 $(px+(1-p))^n \text{ mod } x^{c_1}-1$ 即可。

拉格朗日反演

有两个多项式 $F(x)$ 和 $G(x)$，均满足常数项为 $0$ 且 $1$ 次项不为 $0$ ，且两者为复合逆 $G(F(x))=1$。由拉格朗日反演可以得到：

$$[x^n]F(x)=\frac{1}{n}[x^{-1}]\frac{1}{G(x)}$$

在算法竞赛中我们更关心这个形式：

$$[x^n]F(x)=\frac{1}{n}[x^{n-1}](\frac{x}{G(x)})^{n}$$

BZOJ3684 大朋友和多叉树

题意

对于一棵带有正整数点权的有根多叉树，如果它满足这样的性质，我们的大朋友就会将其称作神犇的：点权为 $1$ 的结点是叶子结点；对于任一点权大于 $1$ 的结点 $u$ ，$u$ 的孩子数目 $\text{deg}_u$ 属于集合 $D$ ，且 $u$ 的点权等于这些孩子结点的点权之和。
给出一个整数 $s$ ，你能求出根节点权值为 $s$ 的神犇多叉树的个数。

题解

记答案的 OGF 为 $T(x)$ ，可以写出：

$$\begin{aligned} T(x)=x+\prod_{d\in D} T^d(x) T(x)-\prod_{d\in D} T^d(x)=x\\ \end{aligned}$$

令 $F(x)=T(x)$ , $G(x) = x-\prod_{d\in D} x^d$ 可以凑出拉格朗日反演的形式：

$$G(F(x))=x$$

那么可以通过：

$$[x^n]F(x)=\frac{1}{n}[x^{n-1}](\frac{x}{G(x)})^{n}$$

求逆加倍增在 $O(n\log n)$ 的复杂度求出 $[x^n]F(x)$。