伯努利数学习笔记

下 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……后来重写了一遍就过了。

复制

#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;
}