[UOJ 50][UR #3]链式反应
给你一个集合\(S\)（保证其中元素都为小于\(n\)的自然数），定义一棵合法的树为一棵编号满足堆的性质，且非叶子节点都一定有两个可能非叶子的儿子，同时有\(c(c\in S)\)个一定为叶子的儿子。对于所有\(i = 1\ldots n\)，求有多少大小为\(i\)的形态不同的合法的树。
\(n\leq 200000\)。
毒瘤，毒瘤.jpg
先考虑DP，定义状态\(f(i)\)表示有\(i\)个点的合法树的数量，那么有：
\[
\begin{aligned}
f(n) &= \frac{1}{2}\sum_{j}\sum_{k}[(n - 1 - j - k)\in S]\binom{n - 1}{j}\binom{n - 1 - j}{k}f(j)f(k)\\
&= \frac{1}{2}\sum_{j}\sum_{k}[(n - 1 - j - k)\in S]\frac{(n - 1)!}{(n - 1 - j - k)!}\frac{f(j)}{j!}\frac{f(k)}{k!}
\end{aligned}
\]
那么定义\(g(n) = \sum_{i + j = n} \frac{f(i)f(j)}{i!j!}\)，那么我们可以直接枚举\(j + k\)，所以转移方程变成了：
\[
\begin{aligned}
f(n) &= \frac{1}{2}\sum_{s}[(n - 1 - s)\in S]\frac{(n - 1)!}{(n - 1 - s)!}g(s)
\end{aligned}
\]
然后就能\(O(n^2)\)做啦！虽然会T飞……
然后这里面出来了阶乘和组合数，考虑上指数型生成函数（下称EGF）！定义\(F(x)\)为\(f(i)\)的EGF，那么\(g(i)\)的EGF就是\(F^2(x)\)，我们再定义\(C(x) = \sum_{i = 0}^{+\infty} [i\in S]\frac{1}{i!}\)，然后我们得到了这样一个关系：
\[\frac{1}{2}C(x)F^2(x) + 1 = \frac{\mathrm{d}F(x)}{\mathrm{d}x}\]
然后你会发现这个微分方程搞不动了……（就是怎么查解都查不出来）
那么上多项式牛迭？我们还是考虑从膜低次的情况（假设为\(F_0(x)\)）推到高次情况（假设为\(F(x)\)，此时模式为\(x^p\)），顺便设\(h(t) = \frac{1}{2}C(x)t + 1\)，那么很显然有\(\frac{\mathrm{d}F(x)}{\mathrm{d}x}\equiv h(F(x))\pmod{x^p}\)，那这怎么搞？
我们还是使用传统艺能，将\(h\)在\(F_0(x)\)处泰勒展开，那么得到：
\[
\begin{aligned}
h(F(x))&\equiv h(F_0(x)) + \frac{h'(F_0(x))}{1!}(F(x) - F_0(x))^1 +\ldots\pmod{x^p}\\
&\equiv h(F_0(x)) + h'(F_0(x))F(x) - h'(F_0(x))F_0(x)\pmod{x^p}\\
\frac{\mathrm{d}F(x)}{\mathrm{d}x}&\equiv h(F_0(x)) + h'(F_0(x))F(x) - h'(F_0(x))F_0(x)\pmod{x^p}
\end{aligned}
\]
然后这泥萌还是没法做啊……
我们再挣扎一下，设\(v(x) = e^{-\int h'(F_0(x))\mathrm{d}x}\)，那么我们将上面的等式两边乘\(v(x)\)，得：
\[
\begin{aligned}
\frac{\mathrm{d}F(x)}{\mathrm{d}x}v(x)&\equiv v(x)(h(F_0(x)) + h'(F_0(x))F(x) - h'(F_0(x))F_0(x))\pmod{x^p}\\
\frac{\mathrm{d}F(x)}{\mathrm{d}x}v(x)-h'(F_0(x))F(x)v(x)&\equiv v(x)(h(F_0(x)) - h'(F_0(x))F_0(x))\pmod{x^p}\\
\frac{\mathrm{d}F(x)}{\mathrm{d}x}v(x)+F(x)\frac{\mathrm{d}v(x)}{\mathrm{d}x}&\equiv v(x)(h(F_0(x)) - h'(F_0(x))F_0(x))\pmod{x^p}\\
\frac{\mathrm{d}(F(x)v(x))}{\mathrm{d}x}&\equiv v(x)(h(F_0(x)) - h'(F_0(x))F_0(x))\pmod{x^p}\\
F(x)v(x)&\equiv \int v(x)(h(F_0(x)) - h'(F_0(x))F_0(x))\mathrm{d}x\pmod{x^p}\\
F(x)&\equiv \frac{1}{v(x)}\int v(x)(h(F_0(x)) - h'(F_0(x))F_0(x))\mathrm{d}x\pmod{x^p}
\end{aligned}
\]
然后这TM就可以倍增了……至于多项式exp、积分、求逆都是传统艺能了……顺便注意下常数……
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cassert>
#include <algorithm>
#include <functional>
#include <utility>
using ll = long long;
constexpr ll ha = 998244353LL;
inline ll pow_mod(ll a, ll b) {
ll ans = 1, res = a;
while(b) {
if(1LL & b) ans = ans * res % ha;
res = res * res % ha; b >>= 1;
}
return ans;
}
inline ll inv(ll x) {
return pow_mod(x, ha - 2LL);
}
const int maxn = 524300;
ll fac[maxn], ifac[maxn], invp[maxn];
inline void process() {
invp[1] = 1;
for(int i = 2; i < maxn; i ++) {
invp[i] = (ha - (ha / (ll)i));
invp[i] = invp[i] * invp[ha % (ll)i] % ha;
}
fac[0] = 1;
for(int i = 1; i < maxn; i ++) {
fac[i] = fac[i - 1] * (ll)i % ha;
}
ifac[maxn - 1] = inv(fac[maxn - 1]);
for(int i = maxn - 2; i >= 0; i --) {
ifac[i] = (ll(i + 1)) * ifac[i + 1] % ha;
}
}
inline void NTT(ll *A, int bi, bool flag = false) {
int n = 1 << bi;
for(int i = 1, j = 0; i < n; i ++) {
int k = n;
do j ^= (k >>= 1); while ((j & k) == 0);
if(i < j) std::swap(A[j], A[i]);
}
for(int L = 1; L < n; L <<= 1) {
const bool bb = (L < 4);
ll c = (ha - 1LL) / (ll(L << 1));
if(flag) c = ha - 1LL - c;
ll xi_n = pow_mod(3LL, c);
ll *B = A + L;
for(int i = 0; i < n; i += (L << 1)) {
ll xi = 1;
if(true) {
for(int j = i; j < i + L; j ++) {
ll x = A[j], y = B[j] * xi % ha;
A[j] = (x + y) % ha;
B[j] = (x - y + ha) % ha;
xi = xi * xi_n % ha;
}
} else {
for(int j = i; j < i + L; j += 4) {
ll x = A[j], y = B[j] * xi % ha;
A[j] = (x + y); if(A[j] >= ha) A[j] -= ha;
B[j] = (x - y + ha); if(B[j] >= ha) B[j] -= ha;
xi = xi * xi_n % ha;
x = A[j + 1], y = B[j + 1] * xi % ha;
A[j + 1] = (x + y); if(A[j + 1] >= ha) A[j + 1] -= ha;
B[j + 1] = (x - y + ha); if(B[j + 1] >= ha) B[j + 1] -= ha;
xi = xi * xi_n % ha;
x = A[j + 2], y = B[j + 2] * xi % ha;
A[j + 2] = (x + y); if(A[j + 2] >= ha) A[j + 2] -= ha;
B[j + 2] = (x - y + ha); if(B[j + 2] >= ha) B[j + 2] -= ha;
xi = xi * xi_n % ha;
x = A[j + 3], y = B[j + 3] * xi % ha;
A[j + 3] = (x + y); if(A[j + 3] >= ha) A[j + 3] -= ha;
B[j + 3] = (x - y + ha); if(B[j + 3] >= ha) B[j + 3] -= ha;
xi = xi * xi_n % ha;
}
}
}
}
if(flag) {
ll inv_n = invp[n];
for(int i = 0; i < n; i ++) {
A[i] = A[i] * inv_n % ha;
}
}
}
void poly_inv(int mod, ll *B, ll *BB) {
if(mod == 1) {
BB[0] = inv(B[0]);
} else {
poly_inv((mod + 1) >> 1, B, BB);
int bi = 0, sz = 1;
while(sz <= ((mod * 2) - 1)) {
bi ++; sz <<= 1;
}
static ll tmp[maxn];
std::copy(B, B + mod, tmp);
std::fill(tmp + mod, tmp + sz, 0LL);
NTT(tmp, bi); NTT(BB, bi);
for(int i = 0; i < sz; i ++) {
tmp[i] = (tmp[i] * BB[i]) % ha;
tmp[i] = (tmp[i] * (ha - 1LL)) % ha;
tmp[i] = (tmp[i] + 2LL) % ha;
tmp[i] = (tmp[i] * BB[i]) % ha;
}
NTT(tmp, bi, true);
std::copy(tmp, tmp + mod, BB);
std::fill(BB + mod, BB + sz, 0LL);
}
}
inline void poly_ln(int mod, ll *B, ll *BB) {
int bi = 0, sz = 1;
while(sz <= ((mod * 2) - 1)) {
bi ++; sz <<= 1;
}
static ll tmp[maxn]; memset(tmp, 0, sizeof(ll) * sz);
poly_inv(mod, B, BB);
for(int i = 0; i < mod - 1; i ++) tmp[i] = B[i + 1] * (ll(i + 1)) % ha;
NTT(tmp, bi); NTT(BB, bi);
for(int i = 0; i < sz; i ++) BB[i] = BB[i] * tmp[i] % ha;
NTT(BB, bi, true);
for(int i = mod - 1; i > 0; i --) BB[i] = BB[i - 1] * invp[i] % ha;
BB[0] = 0; std::fill(BB + mod, BB + sz, 0LL);
// free(tmp);
}
void poly_exp(int mod, ll *B, ll *BB) {
if(mod == 1) {
BB[0] = 1;
} else {
poly_exp((mod + 1) >> 1, B, BB);
int bi = 0, sz = 1;
while(sz <= ((mod * 2) - 1)) {
bi ++; sz <<= 1;
}
std::fill(BB + mod, BB + sz, 0LL);
static ll tmp[maxn]; memset(tmp, 0, sizeof(ll) * sz);
poly_ln(mod, BB, tmp);
for(int i = 0; i < mod; i ++) {
tmp[i] = (B[i] - tmp[i] + ha);
if(tmp[i] >= ha) tmp[i] -= ha;
}
tmp[0] ++;
NTT(tmp, bi); NTT(BB, bi);
for(int i = 0; i < sz; i ++) BB[i] = tmp[i] * BB[i] % ha;
NTT(BB, bi, true);
std::fill(BB + mod, BB + sz, 0LL);
// free(tmp);
}
}
void poly_equation(int mod, ll *B, ll *BB) {
if(mod == 1) {
BB[0] = 0;
} else {
poly_equation((mod + 1) >> 1, B, BB);
int bi = 0, sz = 1;
while(sz <= (mod * 2 - 1)) {
bi ++; sz <<= 1;
}
static ll t1[maxn], rv[maxn], v[maxn];
std::copy(B, B + mod, t1); memset(t1 + mod, 0, (sz - mod) * sizeof(ll));
NTT(t1, bi);
// std::fill(BB + mod, BB + sz, 0LL);
NTT(BB, bi);
for(int i = 0; i < sz; i ++) rv[i] = (t1[i] * BB[i] % ha);
NTT(rv, bi, true);
for(int i = mod - 1; i > 0; i --) rv[i] = (ha - rv[i - 1] * invp[i] % ha);
rv[0] = 0;
poly_exp(mod, rv, v);
static ll rn[maxn]; memset(rn, 0, sizeof(ll) * sz);
// for(int i = 0; i < sz; i ++) rn[i] = BB[i] * BB[i] % ha;
// NTT(rn, bi, true); std::fill(rn + mod, rn + sz, 0LL);
// NTT(rn, bi);
static const ll ii_2 = ha - inv(2);
for(int i = 0; i < sz; i ++) {
rn[i] = (1LL + (ii_2 * (BB[i] * BB[i] % ha) % ha) * t1[i]) % ha;
}
NTT(rn, bi, true); std::fill(rn + mod, rn + sz, 0LL);
memset(rv, 0, sizeof(ll) * sz);
poly_inv(mod, v, rv); NTT(v, bi); NTT(rn, bi);
for(int i = 0; i < sz; i ++) rn[i] = rn[i] * v[i] % ha;
NTT(rn, bi, true); std::fill(rn + mod, rn + sz, 0LL);
for(int i = mod - 1; i > 0; i --) rn[i] = rn[i - 1] * invp[i] % ha;
rn[0] = 0;
NTT(rn, bi); NTT(rv, bi);
for(int i = 0; i < sz; i ++) rn[i] = rn[i] * rv[i] % ha;
NTT(rn, bi, true);
std::copy(rn, rn + mod, BB); std::fill(BB + mod, BB + sz, 0LL);
// free(rn); free(v); free(rv); free(t1);
}
}
char S[maxn];
ll A[maxn], C[maxn];
int main() {
#ifdef LOCAL
freopen("gou.in", "r", stdin);
freopen("outp", "w", stdout);
srand(238238);
for(int i = 0; i < 200000; i ++) {
if(rand() & 1) S[i] = '1';
}
#endif
process();
int n; scanf("%d%s", &n, S);
for(int i = 0; i < n; i ++) {
if(S[i] == '1') {
C[i] = ifac[i];
}
}
poly_equation(n + 1, C, A);
for(int i = 1; i <= n; i ++) {
printf("%lld\n", A[i] * fac[i] % ha);
}
return 0;
}
评论 (0)
