伯努利数学习笔记

发表于 | 更新于 | 分类于 | 阅读次数：

下 Summer Pockets REFLECTION BLUE 的时候太无聊，于是就来学这个东西了。

伯努利数(Bernoulli number)是由雅各布·伯努利的名字命名的。

定义

我们定义伯努利数 $B_n$ 是满足如下递归式 $\sum_{i=0}^{n}\binom{n + 1}{i}B_i = [n = 0]$ 的数列。
定义自然数幂和函数 $S(n, k)=\sum\limits_{i=0}^{n - 1}i^k$。

EGF

首先求出它的 EGF $B(z)$。定义式两边同加 $B_{n + 1}$ 得 $\sum\limits_{i = 0}^{n + 1} \binom{n + 1}{i} B_i = [n = 0] + B_{n + 1}$。

处理一下就是 $\sum\limits_{i = 0}^{n} \binom{n}{i} B_i = [n = 1] + B_{n}$。

所以 $B(z)e^z = B(z) + z$，即 $B(z) = \dfrac{z}{e^z - 1}$。

与自然数幂的转化

$$S(n, k) = \sum_{i = 0}^{n - 1} i^k = \dfrac{1}{k + 1}\sum_{i = 0}^{k}\binom{k + 1}{i} B_i n^{k - i + 1}$$

这玩意大概有两种常见证法。

利用归纳法证明

这个证明方法来自 Concrete Mathematics 6.5 BERNOULLI NUMBER。

$$\begin{aligned}
S(n, k)+n^{k + 1} &= \sum_{i = 0}^{n - 1}(i + 1)^{k + 1}\\
&= \sum_{i = 0}^{n - 1}\sum_{j=0}^{k + 1}\binom{k + 1}{j}k^j\\
&= \sum_{j = 0}^{k + 1}\binom{k + 1}{j}S(n, j)
\end{aligned}$$

两边同时减去 $S(n, k)$ 得

$$n^{k + 1} = \sum_{j = 0}^{k} \binom{k + 1}{j} S(n, j)$$

设 $\hat{S}(n, k) = \sum_{i = 0}^{n - 1} i^k = \dfrac{1}{k + 1}\sum_{i = 0}^{k}\binom{k + 1}{i} B_i n^{k - i + 1}$，且对于 $j\in [0, k)$ 均有 $S(n, j) = \hat{S}(n, j)$ 成立。则需证明 $S(n, k) = \hat{S}(n, k)$。

则此时 $n^{k + 1} = \sum_{j = 0}^{k - 1} \binom{k + 1}{j} \hat{S}(n, j) + S(n, k)$，故转化为证 $n^{k + 1} = \sum_{j = 0}^{k} \binom{k + 1}{j}\hat{S}(n, j)$。

$$\begin{aligned}
&\sum_{j = 0}^{k} \binom{k + 1}{j}\hat{S}(n, j)\\
&= \sum_{j = 0}^{k} \binom{k + 1}{j}\dfrac{1}{j + 1}\sum_{i = 0}^{j}\binom{j + 1}{i} B_i n^{j - i + 1}\\
&= \sum_{j = 0}^{k} \binom{k + 1}{j}\dfrac{1}{j + 1}\sum_{i = 0}^{j}\binom{j + 1}{i + 1} B_{j - i} n^{i + 1}\\
&= \sum_{j = 0}^{k} \binom{k + 1}{j}\sum_{i = 0}^{j}\binom{j}{i} B_{j - i} \dfrac{n^{i + 1}}{i + 1}\\
&= \sum_{i = 0}^{k} \dfrac{n^{i + 1}}{i + 1}\sum_{j = i}^{k}\binom{k + 1}{j}\binom{j}{i} B_{j - i}\\
&= \sum_{i = 0}^{k} \dfrac{n^{i + 1}}{i + 1}\binom{k + 1}{i}\sum_{j = i}^{k}\binom{k + 1 - i}{j - i} B_{j - i}\\
&= \sum_{i = 0}^{k} \dfrac{n^{i + 1}}{i + 1}\binom{k + 1}{i}\sum_{j = 0}^{k - i}\binom{k + 1 - i}{j} B_{j}\\
&= \sum_{i = 0}^{k} \dfrac{n^{i + 1}}{i + 1}\binom{k + 1}{i}[i = k]\\
&= \dfrac{n^{k + 1}}{k + 1}\binom{k + 1}{k}\\
&= n^{k + 1}
\end{aligned}$$

其中倒数第二步用的是定义式。

利用指数生成函数证明

设 $F_n(z) = \sum\limits_{k\ge 0}\dfrac{S(n, k)}{k!}z^k$。

$$\begin{aligned}
F_n(z) &= \sum\limits_{k\ge 0}\dfrac{S(n, k)}{k!}z^k\\
&= \sum_{i=0}^{n-1}\sum_{k\ge 0}\dfrac{i^kz^k}{k!}\\
&= \sum_{i=0}^{n-1}e^{iz}\\
&= \dfrac{e^{nz} - 1}{e^z - 1}\\
&= \dfrac{z}{e^z - 1}\cdot\dfrac{e^{nz} - 1}{z}\\
&= B(z)\cdot\dfrac{e^{nz} - 1}{z}\\
&= \left(\sum_{i\ge 0}\dfrac{B_i}{i!} \right)\left(\sum_{i\ge 0}\dfrac{n^{i+1} z^{i}}{(i+1)!}\right)
\end{aligned}$$

然后取一项就行了。

$$\begin{aligned}
S(n, k) &= k![z^k]F_n(z)\\
&= m!\sum_{i=0}^{k}\dfrac{B_i}{i!}\cdot\dfrac{n^{k-i+1}}{(k-i+1)!}\\
&= \dfrac{1}{k+1}\sum_{i=0}^{k}\binom{k+1}{i}B_in^{k-i+1}
\end{aligned}$$

例题

未特别说明时，模数均为998244353。

「Luogu 3711」仓鼠的数学题

Luogu 3711

给数列 ${a_n}$，求 $Ans(z) = \sum\limits_{k = 0}^{n}(S(z, k) + z^k)a_k$。$n\le 2.5\times 10^5$。

$\sum\limits_{k = 0}^n a_k z^k$ 是个常量可以先扔出来。

$$\begin{aligned}
&Ans(z) -\sum\limits_{k = 0}^n a_k z^k\\
&= \sum_{k = 0}^{n}S(z, k)a_k\\
&= \sum_{k = 0}^{n}a_k\dfrac{1}{k+1}\sum_{i=0}^{k}\binom{k + 1}{i}B_i z^{k - i + 1}\\
&= \sum_{k = 0}^{n}a_k k!\sum_{i = 0}^{k} \dfrac{B_i}{i!}\cdot\dfrac{z^{k + 1 - i}}{(k + 1 - i)!}\\
&= \sum_{k = 0}^{n}a_k k!\sum_{i = 1}^{k + 1} \dfrac{B_{k + 1 - i}}{(k + 1 - i)!}\cdot\dfrac{z^i}{i!}\\
&= \sum_{i = 1}^{n + 1}\dfrac{z^i}{i!}\sum_{k = i - 1}^{n}(a_k k!)\dfrac{B_{k + 1 - i}}{(k + 1 - i)!}
\end{aligned}$$

后面那坨是个差卷积，于是翻转一下多项式做完了。

「HDU 6340」Delightful Formulas

HDU 6340

给一个大数 $n = \prod\limits_{o = 0}^{m - 1} p_o^{\alpha_o}$ 的分解和正整数 $K$，求 $\sum\limits_{i=1}^{n}[\gcd(i, n) = 1]S(i, K)$。
$m\le 20, K\le 10^5, p_o, \alpha_o\le 10^9$。

首先来个莫比乌斯反演，$Ans = \sum\limits_{d|n} \mu(d)\sum\limits_{i=1}^{n/d}S(i\cdot d, K)$。设 $a_i = \dfrac{1}{k + 1}\binom{k + 1}{j}B_{k + 1 - j}$。

$$\begin{aligned}
Ans &= \sum_{d|n}\mu(d)\sum_{j=1}^{K+1}a_j d^j \sum_{i = 1}^{n/d} i^j\\
&= \sum_{d|n}\mu(d)\sum_{j=1}^{K+1}a_j d^j \sum_{i=1}^{j+1}(\dfrac{n}{d})^i\dfrac{1}{j + 1}\binom{j + 1}{i}B_{j + 1 - i}
\end{aligned}$$

到这里，枚举 $c = j - i + 1\in[0, K + 1]$：

$$\begin{aligned}
Ans &= \sum_{d|n}\mu(d)\sum_{c=-}1^{K}d^c\sum_{j=1}^{K+1}a_j n^{j-c} \dfrac{1}{j + 1}\binom{j+1}{c+1}B_{c+1}\\
&= \sum_{c=-1}^{K}B_{c+1}\sum_{j=1}^{K+1}\dfrac{a_j}{j+1}\binom{j+1}{c+1}n^{j-c} \sum_{d|n}d^c\mu(d)\\
&= \sum_{c=0}^{K+1}B_{c}\sum_{j=1}^{K+1}\dfrac{a_j}{j+1}\binom{j+1}{c}n^{j-c+1} \sum_{d|n}d^{c-1}\mu(d)
\end{aligned}$$

发现 $\sum\limits_{d|n}d^c\mu(d)=\prod\limits_{o=0}(1 - p_o)$，所以考虑 $F_c = \sum\limits_{j=1}^{K+1}\dfrac{a_j}{j+1}\binom{j+1}{c}n^{j-c+1}$，我们要对于 $c\in[0, K + 1]$ 全求出来。

$$\begin{aligned}
F_c &= \sum_{j=1}^{K+1}\dfrac{a_j}{j+1}\binom{j+1}{c}n^{j-c+1}\\
&= \dfrac{1}{c!}\sum_{j=1}^{K+1}a_j j!\dfrac{n^{j - c + 1}}{(j - c + 1)!}
\end{aligned}$$

又是差卷积，稍微算算发现 $K-j+c\in[-K-1, K]$，所以设 $b_i=\dfrac{n^{K - i + 2}}{(K - i + 2)!}$。
就变成求 $F_c = \sum\limits_{j=1}^{K+1}a_j j! b_{K + 1-j+c}$。（这个 $a_i$ 其实拆开形式比较好）

于是做完了。记得乘上一堆麻烦的系数。

不知道为什么板子一直 WA……后来重写了一遍就过了。

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
167
168
169
170
171
172
173
174
175
176
177
178
179
#include <bits/stdc++.h>
namespace my_std {
using namespace std;
#define reg register
#define Rint register int
#define FOR(i, a, b) for (register int i = (a), ed_##i = (b); i <= ed_##i; ++i)
#define ROF(i, a, b) for (register int i = (a), ed_##i = (b); i >= ed_##i; --i)
typedef long long i64;
typedef unsigned long long u64;
#define Templ(T) template <typename T>
inline int read() {
    reg int ans = 0, f = 1;
    reg char c = getchar();
    while (!isdigit(c))
        f ^= (c == '-'), c = getchar();
    for (; isdigit(c); c = getchar())
        ans = (ans << 1) + (ans << 3) + (c ^ 48);
    return f ? ans : -ans;
}
Templ(_Tp) inline int chkmin(_Tp &x, _Tp y) { return x > y ? x = y, 1 : 0; }
Templ(_Tp) inline int chkmax(_Tp &x, _Tp y) { return x < y ? x = y, 1 : 0; }
#define using_mod
// const int mod = 998244353;
#define mod 998244353
#ifdef using_mod
inline void inc(int &x, const int &y) {
    x += y;
    if (x >= mod)
        x -= mod;
}
inline void dec(int &x, const int &y) {
    x -= y;
    if (x < 0)
        x += mod;
}
inline int qmo(reg const int &x) { return x + ((x >> 31) & mod); }
inline int ksm(int x, i64 y) {
    int res = 1;
    for (; y; y >>= 1, x = (i64)x * x % mod)
        if (y & 1) res = (i64)res * x % mod;
    return res;
}
#endif
} // namespace my_std
using namespace my_std;

#define swap(x, y) (x ^= y ^= x ^= y)
const int N = 550010;

int LMT = 1;
int rev[N], omg[N], inv[N];
int fac[N], ifac[N];
int l2g[N << 1];
inline void init(const int &n) {
    fac[0] = fac[1] = ifac[0] = ifac[1] = inv[1] = 1;
    FOR(i, 2, n) {
        inv[i] = (i64)(mod - mod / i) * inv[mod % i] % mod;
        fac[i] = (i64)fac[i - 1] * i % mod;
        ifac[i] = (i64)ifac[i - 1] * inv[i] % mod;
    }
    l2g[1] = 0;
    FOR(i, 2, n << 1)
    l2g[i] = l2g[i >> 1] + 1;
}

inline void poly_init(const int &n){
    Rint l = 0;
    while(LMT <= n) LMT <<= 1, ++l;
    FOR(i, 1, LMT - 1) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
    reg const int t = ksm(3, (mod - 1) >> l);
    omg[LMT >> 1] = 1;
    FOR(i, (LMT >> 1) + 1, LMT - 1) omg[i] = (i64)omg[i - 1] * t % mod;
    ROF(i, (LMT >> 1) - 1, 1) omg[i] = omg[i << 1];
    LMT = l;
}

inline int get_len(const int &n){
    return 1 << (l2g[n] + 1);
}

inline void DFT(int *a, const int &n){
    static u64 tmp[N];
    reg const int fix = LMT - l2g[n];
    Rint t;
    FOR(i, 0, n - 1) tmp[i] = a[rev[i] >> fix];
    for(Rint i = 1; i < n; i <<= 1){
        for(Rint j = 0; j < n; j += i << 1){
            FOR(k, j, j + i - 1){
                t = tmp[i + k] * omg[i + k - j] % mod;
                tmp[i + k] = tmp[k] + mod - t;
                tmp[k] += t;
            }
        }
    }
    FOR(i, 0, n - 1) a[i] = tmp[i] % mod;
}
inline void IDFT(int *a, const int &n){
    reverse(a + 1, a + n);
    DFT(a, n);
    reg const int bk = mod - (mod - 1) / n;
    FOR(i, 0, n - 1) a[i] = (i64)bk * a[i] % mod;
}

inline void poly_mul(int *a, int *b, int *c, const int &deg){
    static int tmp1[N], tmp2[N];
    reg const int len = get_len(deg);
    memcpy(tmp1, a, sizeof(int) * len), memcpy(tmp2, b, sizeof(int) * len);
    DFT(tmp1, len), DFT(tmp2, len);
    FOR(i, 0, len - 1) c[i] = (i64)tmp1[i] * tmp2[i] % mod;
    IDFT(c, len);
    memset(c + deg, 0, sizeof(int) * (len - deg));
}
inline void poly_inv(int *a, int *b, const int &deg){
    static int tmp[N];
    if(deg == 1){
        *b = ksm(*a, mod - 2);
        return;
    }
    poly_inv(a, b, (deg + 1) >> 1);
    reg const int len = get_len(deg << 1);
    memcpy(tmp, a, sizeof(int) * deg);
    memset(tmp + deg, 0, sizeof(int) * (len - deg));
    DFT(b, len), DFT(tmp, len);
    FOR(i, 0, len - 1){
        b[i] = (i64)qmo(2ll - (i64)b[i] * tmp[i] % mod) * b[i] % mod;
    }
    IDFT(b, len);
    memset(b + deg, 0, sizeof(int) * (len - deg));
}

#define M (1 << 17)
int bnl[N], A[N], G[N];
int mc, n, K;
int pr[20], aph[20], pt[20];

inline void work(){
    K = read(), mc = read(), n = 1;
    FOR(o, 0, mc - 1){
        pr[o] = read(), aph[o] = read();
        n = n * (i64)ksm(pr[o], aph[o]) % mod;
    }
    reg const int len = get_len((K + 5) * 2);
    memset(A, 0, sizeof(int) * len);
    memset(G, 0, sizeof(int) * len);
    FOR(i, 0, K){
        A[K + 1 - i] = (i64)ifac[i] * bnl[i] % mod;
    }
    for(Rint i = 1, j = 1; i <= K + 2; ++i){
        j = j * (i64)n % mod;
        G[K + 2 - i] = j * (i64)ifac[i] % mod;
    }
    poly_mul(A, G, G, len - 1);
    Rint ans(0);
    FOR(o, 0, mc - 1) pt[o] = ksm(pr[o], mod - 2);
    FOR(c, 0, K + 1){
        Rint res = G[c + K + 1] * (i64)bnl[c] % mod * ifac[c] % mod;
        FOR(o, 0, mc - 1){
            res = (i64)res * (1ll - pt[o] + mod) % mod;
            pt[o] = pt[o] * (i64)pr[o] % mod;
        }
        ans = qmo(ans + res - mod);
    }
    ans = (i64)fac[K] * ans % mod;
    printf("%d\n", ans);
}

int main() {
    init(M << 1);
    poly_init(M << 1);
    memcpy(G, ifac + 1, sizeof(int) * M);
    poly_inv(G, bnl, M);
    FOR(i, 1, M - 1) bnl[i] = (i64)bnl[i] * fac[i] % mod;
    bnl[1] = mod - bnl[1];

    int esac(read());
    while (esac--)
        work();
    return 0;
}