[BZOJ3684]大朋友和多叉树

[BZOJ3684]大朋友和多叉树

试题描述

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

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

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

输入

第一行有 (2) 个整数 (s,m)。

第二行有 (m) 个互异的整数，(d[1],d[2], cdots ,d[m])，为集合 (D) 中的元素。

输出

输出一行仅一个整数，表示答案模 (950009857) 的值。

输入示例

4 2
2 3

输出示例

10

样例解释

数据规模及约定

(1 le m < s le 10^5, 2 le d[i] le s)，有 (3) 组小数据和 (3) 组大数据。

题解

我们还从暴力 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)，那么就能得到这样一个方程：

egin{equation}
F(x) = x + sum_{j in D} F(x)^j
end{equation}

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

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

(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)} ight) ^n ]

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

可能对上面的公式你还不明白 ([x^{-1}]) 是什么操作。注意到 (g(x)) 的常数项为 (0)，所以它的逆元是不存在的，即没有任何整式能够表示 (frac{1}{g(x)})，这时候就需要扩充一下这个域，引入分式域(^{[3]})。在这里所有的多项式能够被表示成 (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)} ight) ^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) 这个多项式，问题就解决了。

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
using namespace std;
#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)

int read() {
	int x = 0, f = 1; char c = getchar();
	while(!isdigit(c)){ if(c == '-') f = -1; c = getchar(); }
	while(isdigit(c)){ x = x * 10 + c - '0'; c = getchar(); }
	return x * f;
}

#define maxn 262144
#define MOD 950009857
#define Groot 7
#define LL long long

int Pow(int a, int b) {
	int ans = 1, t = a;
	while(b) {
		if(b & 1) ans = (LL)ans * t % MOD;
		t = (LL)t * t % MOD; b >>= 1;
	}
	return ans;
}

int inv[maxn];
void init() {
	inv[1] = 1;
	rep(i, 2, maxn - 1) inv[i] = (LL)(MOD - MOD / i) * inv[MOD%i] % MOD;
	return ;
}

int brev[maxn];
void FFT(int *a, int len, int tp) {
	int n = 1 << len;
	rep(i, 0, n - 1) if(i < brev[i]) swap(a[i], a[brev[i]]);
	rep(i, 1, len) {
		int wn = Pow(Groot, MOD - 1 >> i);
		if(tp < 0) wn = Pow(wn, MOD - 2);
		for(int j = 0; j < n; j += 1 << i) {
			int w = 1;
			rep(k, 0, (1 << i >> 1) - 1) {
				int la = a[j+k], ra = (LL)w * a[j+k+(1<<i>>1)] % MOD;
				a[j+k] = (la + ra) % MOD;
				a[j+k+(1<<i>>1)] = (la - ra + MOD) % MOD;
				w = (LL)w * wn % MOD;
			}
		}
	}
	if(tp < 0) {
		int invn = Pow(n, MOD - 2);
		rep(i, 0, n - 1) a[i] = (LL)a[i] * invn % MOD;
	}
	return ;
}

void Mul(int *A, int *B, int n, int m, bool recover = 0) {
	int N = 1, len = 0;
	while(N <= n + m) N <<= 1, len++;
	rep(i, 0, N - 1) brev[i] = (brev[i>>1] >> 1) | ((i & 1) << len >> 1);
	rep(i, n + 1, N - 1) A[i] = 0;
	rep(i, m + 1, N - 1) B[i] = 0;
	FFT(A, len, 1); FFT(B, len, 1);
	rep(i, 0, N - 1) A[i] = (LL)A[i] * B[i] % MOD;
	FFT(A, len, -1); if(recover) FFT(B, len, -1);
	return ;
}

int tmp[maxn];
void inverse(int *f, int *g, int n) {
	if(n == 1) return (void)(f[0] = Pow(g[0], MOD - 2));
	inverse(f, g, n + 1 >> 1);
	rep(i, 0, n - 1) tmp[i] = g[i];
	int N = 1, len = 0;
	while(N < (n << 1)) N <<= 1, len++;
	rep(i, 0, N - 1) brev[i] = (brev[i>>1] >> 1) | ((i & 1) << len >> 1);
	rep(i, n, N - 1) tmp[i] = 0;
	rep(i, n + 1 >> 1, N - 1) f[i] = 0;
	FFT(f, len, 1); FFT(tmp, len, 1);
	rep(i, 0, N - 1) f[i] = ((LL)f[i] * (2ll - (LL)tmp[i] * f[i] % MOD) % MOD + MOD) % MOD;
	FFT(f, len, -1);
	return ;
}

int deri[maxn], pinv[maxn];
void logarithm(int *f, int *g, int n) { // g[0] must be 1
	rep(i, 1, n - 1) deri[i-1] = (LL)g[i] * i % MOD;
	inverse(pinv, g, n);
	Mul(deri, pinv, n - 2, n - 1);
	rep(i, 0, n - 2) f[i+1] = (LL)deri[i] * inv[i+1] % MOD; f[0] = 0;
	return ;
}

int lnf[maxn], tg[maxn];
void exponential(int *f, int *g, int n) { // g[0] must be 0
	if(n == 1) return (void)(f[0] = 1);
	exponential(f, g, n + 1 >> 1);
	rep(i, n + 1 >> 1, n - 1) f[i] = 0;
	logarithm(lnf, f, n);
	rep(i, 0, n - 1) tg[i] = (g[i] - lnf[i] + MOD) % MOD; tg[0]++; if(tg[0] >= MOD) tg[0] -= MOD;
	Mul(f, tg, (n + 1 >> 1) - 1, n - 1);
	return ;
}

int t1[maxn];
void p_pow(int *a, int n, int b) { // a[0] must be 1
	logarithm(t1, a, n);
	rep(i, 0, n - 1) t1[i] = (LL)t1[i] * b % MOD;
	exponential(a, t1, n);
	return ;
}

int n, m, d[maxn], ans[maxn];
int main() {
	init();
	
	n = read(); m = read();
	rep(i, 1, m) d[read()-1] = MOD - 1;
	d[0] = 1;
	p_pow(d, n, n);
	inverse(ans, d, n);
	
	printf("%lld
", (LL)inv[n] * ans[n-1] % MOD);
	
	return 0;
}

([1 sim 3])：拉格朗日反演定理的证明，为什么 (f(g(x)) = g(f(x)) = x) 以及所谓的“分式域”：本文为尽量通俗易懂地阐述此题解法并没有严格地用数学符号，也没有给出一些结论的证明，详见这里。