「算法学习」多项式与生成函数

我 就 是 机 械 毁 灭 者

写在开头:写多项式最重要的事情就是注意清空,不仅仅是低次项或是两倍长度的低次项,就算是有很高次项没有清空,也会发生一些匪夷所思的事情,毕竟 NTT 这种东西没什么道理好讲,如果有数组贴着极限开的习惯可以直接在牛顿迭代的时候 memset,不会影响复杂度,否则在做完每一次 ntt 后最好都把高次项没用的丢掉,以防万一。

普通生成函数 OGF

卡特兰数

卡特兰数的第 $n$ 项表示 $n$ 个节点的二叉树数量(不代标号,区分左右子树)。记卡特兰数的第 $n$ 项为 $f_n$ ,那么 $f(0) = 1$,且对于 $n \ge 1$ 有递推式:

记 $F(x)$ 表示卡特兰数的生成函数,那么我们有:

展开这个式子

令 $G(x) = \sum_{k = 0}^{\infty}\binom{\frac12}{k}x^k$

将 $x = -4$ 带入

无序分拆数

分拆数 $P_n$ 表示把整数 $n$ 分拆成若干个无序的数的和的方案数

我们令 $F$ 为 $P_n$ 的生成函数,则 $F = \prod_{k=1}^{\infty}\frac1{(1-x^k)}$

我们可以把 $\frac{1}{1-x^k}$ 这一项理解成在某一种方案中,$k$ 这个数被选了几次。

首先我们把等式两边同时取对数(化乘为加),得:

考虑我们如何求出 $\ln(1-x^k)$,我们令 $g(x) = ln(1-x)$,则有:

然后我们把 $g’$ 积分回去

至此我们求出了 $\ln(1-x)$,将其带回原式

这个式子我们先枚举 $i$ ,再枚举 $t$ ,可以在调和级数的时间内求出,然后 exp 回去,总复杂度 $O(n \log n)$

【算法】五边形数定理与拆分数 ——litble

【题解】整数划分 O($\sqrt n$)的 Naive 做法

整数拆分 - 维基百科,自由度百科全书

BZOJ 3028 食物

设 $n$ 的答案为 $a_n$

令 $f(x)$ 为 $\{a\}$ 的生成函数

BZOJ 4001 概率论

我们设 $n$ 个节点的所有二叉树的叶子节点总数为 $f_n$ 个,本质不同的二叉树个数为 $g_n$。

则有递推式:

令 $F(x)$ 为 $f_x$ 的生成函数,$G(x)$ 为 $g_x$ 的生成函数,则有:

我们把 $\frac{F(x)}{x}$ 进行不定积分:

说明我们只要求出了 $H(x)$,然后进行求导就可以很方便的求出 $F_x$ ,此处的 $H(x)$ 就是我们上面提到的卡特兰数,于是我们就求出了 $F(x)$

指数生成函数 EGF

对于实数序列 $\{a_n\}$ ,定义它的指数生成函数是无穷级数

引入这个记号仅仅在于简化运算,因为有一下泰勒级数成立:

指数生成函数又一个很好的性质,设有两个EGF $\hat{A}(x), \hat{B}(x)$,则他们的卷积的指数生成函数 $\hat{C}(x) = \hat{A}(x)\hat{B}(x) = \sum_{k = 0}^{\infty}\sum_{i = 0}^{k}\binom{k}{i}A_iB_{k-i}$

poj3734 Blocks

我们直接写出答案的指数生成函数:

BZOJ3456 城市规划

我们设 $f_x$ 为 $x$ 个点的联通图个数,$g_x$ 为 $x$ 个点的图的个数。

显然,$g_x$ 就是每一图中条边可以连或不连,$g_x = 2^{\binom{x}{2}}$

我们还有另一种 $g$ 的表达方式,考虑一个图的联通块个数区间为 $[1, n]$ ,于是我们可以把 $g$ 的 $EGF\; \hat{G}$ 用 $f$ 的 $EGF\;\hat{F}$ 表达:

此处要在指数生成函数的基础上多去掉一个阶乘,原因是每个联通块是等价的,我们先放哪个方案都是一样的。

例如在 $n = 4$ 的时候,我们先放一个 4 的联通块,再放一个 2 的联通块,方案是 $\binom{6}{4}f_4f_2$,与我们先放一个 2 的联通块,再放一个 4 的联通块 方案是 $\binom{6}{2}f_2f_4$ ,是等价的,所以我们在卷积过后还要在去掉一个阶乘。

注意到 $\hat{G}= 1 + \hat{F} + \frac{\hat{F}^2}{2!} + \frac{\hat{F}^3}{3!}\cdots = e^{\hat{F}}$,则 $\hat{F} = \ln \hat{G}$,我们求出 $F$ 的某一项后,记得把阶乘乘回去,此处乘阶乘是因为 $\hat{F}$ 是指数生成函数,与上方的那个阶乘要区分。

51nod 1728 不动点

设 $\hat{F}(x)$ 为答案的 EGF。显然, $\hat{F}(1) = e^x$ 。经过分析,我们发现答案等价于求 $n$ 个点高度不超过 $k$ 的森林的数量。设 $\hat{G}(k)$ 为高度不超过 $k$ 的树的数量的 EGF,可以得到 $\hat{F}(k)=e^{\hat{G}(k)}$ 。考虑如何求 $\hat{G}(k)$ ,可以枚举一个点作为根,那么显然方案就是 $n[x^{n-1}]\hat{F}(k-1)$ ,化简计算一下可得 $\hat{G}(k)=x\hat{F}(k - 1)$ 。因此 $\hat{F}(k)=e^{x\hat{F}(k-1)}$ 。

CF 831E

首先设 $b_i$ 为每个 $a_i$ 最后会被减多少次,通过归纳可以证明,$E(ans) = E(\prod_{i=1}^na_i-\prod_{i=1}^n(a_i-b_i))=\prod_{i=1}^na_i-E(\prod_{i=1}^n(a_i-b_i))$ 。

现在我们想办法把式子的后半段求出,

我们构造 EGF $\hat{F} = \prod_{i=1}^{n}\sum_{j=1}^{\infty}\frac{(a_i - j)}{j!}x^j$

我们构造 OGF $G(x) = \prod_{i}^n(a_i-x) = \sum_{i = 0}^{n}c_ix^i$ ,则答案为

CF 438E 小朋友和二叉树

我们设 $c_i$ 为 $i$ 是否在集合里,$f_i$ 为权值为 $i$ 的树的方案数,可以得到递推式 :

写成生成函数的形式:

此处 1 为 $F_0$ ,即常数项,解方程可得:

在取减号的时候会使式子在 0 处没有意义,所以舍去。

多项式牛顿迭代

泰勒展开

通俗的讲,就是通过 $f(x)$ 与 $f(x)$ 的导数在 $x_0$ 的取值求出 $f(x)$

牛顿迭代 Newton’s Method

求导法则

基本形式

给定函数 $G(z)$ ,求一个多项式 $F(z) \bmod z^n$ ,满足方程

首先 $n = 1$ ,$G(F(z)) \equiv 0\pmod z$ ,这是要单独求出来的。

此处的 $\bmod z^n$ 的意义仅仅表示我们只考虑这个多项式的最低 $n$ 次项,并不是明确的指明是 $z$ 或是 $F(z)$ 什么的,所以不要多想,对于不同的题目边界各不相同,例如开根的边界是 $f[0] = sqrt(G[0])$,求逆的边界是 $F[0] = G[0]^{-1}$ 等等

假设现在我们已经求出了 $G(F_0(z)) \equiv 0 \pmod{z^{\lceil\frac n2\rceil}}$,考虑如何扩展到 $\bmod z^n$。

我们对 $G(F(z))$ 在 $F_0(z)$ 处进行泰勒展开。

我们注意到 $F_0(z)$ 和 $F(z)$ 的最低 $\lceil\frac n2\rceil$ 次项是相同的,所以 $(F(z) - F_0(z))^2$ 的最低 $2\lceil\frac n2\rceil$ 次项必然都是 $0$ ,所以原式 可以写成:

由于 $G(F(z)) \equiv 0 \pmod{z^n}$ ,所以原式又可以写作:

牛顿迭代在多项式问题上的应用

多项式开方

给定 $A(x)$ 求 $B(x)$ ,满足方程 :

移项可得到方 $F^2(x)-A(x) \equiv 0 \pmod{x^n}$

令 $G(F(x)) = F^2(x) - A(x)$ , $G’(F(x)) = 2F(x)$ ,我们已经将这个问题规约到了牛顿迭代的基本形式了。

多项式求逆

给定 $A(x)$ 求 $B(x)$,满足方程:

我们令 $F(x) = B(x)$ , $G(F(x)) = F^{-1}(x) - A(x)$, $G’(F(x)) = -F^{-2}(x)$ ,再次规约到了牛顿迭代的一般形式

多项式的对数和指数函数

我们可以认为多项式的对数是一个多项式和麦克劳林对数的复合,也就是给出一个多项式 $A(x) = \sum_{i\ge1}a_ix^i$

  • 多项式 ln

    对于一个多项式 $A(x) = 1 +\sum_{i\ge1}a_ix^i$,现在要计算其对数 $\ln A(x)$

    考虑直接计算过于麻烦,我们计算其求导后的结果

    多项式求逆我们已经可以在 $O(n\log n)$ 的时间内求出,唯一要做的就是重新积分回去,多项式的积分可以在 $O(n)$ 时间内解决,于是我们可以在 $O(n \log n)$ 的时间内解决多项式 ln

    在多项式的常数项为 $1$ 的时候,我们求导过后再积分,最后常数项通过变换后会变成 $0$ ,但是当常数项不是 $1$ 的时候,问题将变得比较棘手,我们暂不作讨论

  • 多项式 exp

    对于一个多项式 $A(x)=\sum_{i\ge1}a_ix^i$,现在要计算其指数 $e^{A(z)}$。

    考虑我们无法直接计算,利用牛顿迭代。

    令 $F(x) = \exp(A(x)), G(F(x)) = \ln F(X) - A(x), G’(F(x)) = \frac 1 {F(x)}$

    得到递推式

多项式 pow

给定一个多项式 $A(x)$,你现在要计算 $A^k(x)$

有了求指数和求对数的运算后,那么

时间复杂度 $O(n \log n)$ ,可以解决包括多项式开方的计算。

多项式半家桶

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166

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

ll rev[maxn], g = 3;

void ntt(ll *a, ll lim, ll opt)
{
rep(i, 0, lim - 1) if(i < rev[i]) swap(a[i], a[rev[i]]);
for(ll i = 1; i < lim; i <<= 1)
{
ll gn = power(g, (mod - 1) / (i << 1));
if(!opt) gn = power(gn, mod - 2);
for(ll j = 0; j < lim; j += (i << 1))
{
ll t = 1;
for(ll k = 0; k < i; ++ k, (t *= gn) %= mod)
{
ll x = a[j + k], y = a[i + j + k] * t % mod;
a[j + k] = (x + y) % mod;
a[i + j + k] = (x + mod - y) % mod;
}
}
}
if(!opt)
{
ll inv = power(lim, mod - 2);
rep(i, 0, lim - 1) (a[i] *= inv) %= mod;
}
}

void mul(ll *a, ll *b, ll len)
{
ll 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] *= b[i]) %= mod;
ntt(a, lim, 0), ntt(b, lim, 0);
}

void poly_mul(ll *a, ll *b, ll len)
{
ll 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] *= b[i]) %= mod;
ntt(a, lim, 0), ntt(b, lim, 0);
rep(i, len, lim) a[i] = b[i] = 0;
}

ll ta[maxn], tb[maxn];

void dac_mul(ll l, ll r, ll *f, ll *g)
{
if(l == r) return ;
ll 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);
}

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

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

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

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

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

ll inv_ta[maxn], inv_tb[maxn];

void inversion(ll *f, ll *a, ll 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), divide(f, inv_ta, len);
rep(i, 0, len - 1) inv_ta[i] = inv_tb[i] = 0;
}

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

void poly_sqrt(ll *f, ll *a, ll 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;
}

ll der[maxn], ln_inv[maxn];

void poly_ln(ll *f, ll len)
{
rep(i, 0, len - 1) der[i] = ln_inv[i] = f[i];
derivation(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;

}

ll exp_ln[maxn], exp_mul[maxn];

void poly_exp(ll *f, ll *a, ll 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, divide(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;
}

ll prod[maxn];

void poly_pow(ll *f, ll b, ll 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;
}
}

  • 参考文献 :

牛顿迭代法在多项式运算的应用 -Miskcoo’s Space

泰勒级数 -维基百科,自由的百科全书)

多项式部分简介 -OI Wiki

浅谈生成函数之OGF - 知乎

浅谈生成函数之EGF - 知乎

  • Copyright: Copyright is owned by the author. For commercial reprints, please contact the author for authorization. For non-commercial reprints, please indicate the source.
  • Copyrights © 2020-2022 naitir
  • Visitors: | Views:

请我喝杯咖啡吧~

支付宝
微信