伯努利数

定义 $Sk(n)=∑i=0n−1ikS_k(n) = \sum\limits_{i = 0} ^{n - 1} i ^ k$ 。

从二项式出发

$\sum_{i = 0} ^{k} C_{k + 1} ^{i} 0 ^ {i} + 0 ^{k + 1}\\ \vdots\\ (n - 1 + 1) ^{k + 1} = \sum_{i = 0} ^{k} C_{k + 1} ^{i} (n - 1) ^{i} + (n - 1) ^{k + 1}\\ 把次方k + 1的移项，再整体相加，得\\ n ^{k + 1} = \sum_{i = 0} ^{k} C_{k + 1} ^{i} S_i(n)\\ n ^{k + 1} = \sum_{i = 0} ^{k - 1} C_{k + 1} ^{i} S_i(n) + (k + 1) S_k(n)\\ S_k(n) = \frac{1}{k + 1}\left(n ^{k + 1} - \sum_{i = 0} ^{k - 1} C_{k + 1} ^{i} S_i(n) \right)\\$
于是我们有了一个 $O(n ^ 2)$ 递推求 $S_i(n)$ 得算法，

对二项式再进一步考虑：

$^{k} = \sum_{i = 0} ^{k} C_{k} ^{i} 0 ^{i} m ^{k - i}\\ \vdots\\ (n - 1 + m) ^{k} = \sum_{i = 0} ^{k} C_{k} ^{i}(n - 1) ^{i} m ^{k - i}\\ 左右两边同时相加\\ S_k(n + m) - S_k(m) = \sum_{i = 0} ^{k} C_{k} ^{i} S_{i} (n) m ^{k - i}\\ 左右两边对n求导\\ S'_k(n + m) = \sum_{i = 0} ^{k} C_{k} ^{i} S'_i(n) m^{k - i}\\ \\ n = 0，代入\\ S'_k(m) = \sum_{i = 0} ^{k} C_{k} ^{i} S'_i(0) m^{k - i}\\ 两边对m从0->n求积分\\ S_k(n) - S_k(0) = \sum_{i = 0} ^{k} C_{k} ^{i} S'_i(0) \frac{n ^{k - i + 1}}{k - i + 1} - \sum_{i = 0} ^{k} C_{k} ^{i} S'_i(0) \frac{0 ^{k - i + 1}}{k - i + 1}\\ S_k(0) = 0, 对于最右侧i \neq k + 1,所以0 ^{k - i + 1} = 0 \\ S_k(n) = \sum_{i = 0} ^{k} C_{k} ^{i} S'_i(0) \frac{n ^{k - i + 1}}{k - i + 1} = \frac{1}{k + 1} \sum_{i = 0} ^{k} C_{k + 1} ^{i} S'_i(0) n ^{k + i - 1}\\$
其实，这里的 $S_i(0)$ 就是伯努利数，为了方便书写， $B_i = S'_i(0)$ 。
$1代入上式\\ S_k(1) = \frac{1}{k + 1} \sum_{i = 0} ^{k} C_{k + 1} ^{i} B_i\\ 同样的我们把右边的最高项给提出来\\ S_k(1) = \frac{1}{k + 1} \sum_{i = 0} ^{k - 1} C_{k + 1} ^{i} B_i + B_k\\ 则有B_k - S_k(1) = - \frac{1}{k + 1} \sum_{i = 0} ^{k - 1} C_{k + 1} ^{i} B_i\\ 特殊地，考虑k = 0,有S_0(1) = 1, 所以B_0 = 1\\ k \geq 1, S_k(1) = 0, 所以B_k = - \frac{1}{k + 1} \sum_{i = 0} ^{k - 1} C_{k + 1} ^{i} B_i\\$
由此我们可以 $O(k ^ 2)$ 递推出伯努利数出来了。

$F a u l h a b e r$ 公式

$Sk(x)=1k+1∑i=0kCk+1iBixk−i+1S_k(x) = \frac{1}{k + 1} \sum\limits_{i = 0} ^{k}C_{k + 1} ^{i}B_i x ^{k - i + 1}$

一个自然数幂次求和的公式，其实就是我们推导过程中的一步。

#include <bits/stdc++.h>using namespace std;const int N = 3e3 +10, mod = 1e9 + 7;int C[N][N], B[N], inv[N];void init() {for (int i = 0; i < N; i++) {C[i][0] = C[i][i] = 1;for (int j = 1; j < i; j++) {C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % mod;}}inv[1] = 1;for (int i = 2; i < N; i++) {inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;}B[0] = 1;for (int i = 1; i < N - 1; i++) {int cur = 0;for (int j = 0; j < i; j++) {cur = (cur + 1ll * C[i + 1][j] * B[j] % mod) % mod;}cur = 1ll * cur * inv[i + 1] % mod;B[i] = (mod - cur) % mod;}
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);init();int T, k;scanf("%d", &T);while (T--) {long long n, ans = 0, cur = 1;scanf("%lld %d", &n, &k);n++;n %= mod;for (int i = k; i >= 0; i--) {cur = cur * n % mod;ans = (ans + 1ll * C[k + 1][i] * B[i] % mod * cur % mod) % mod;}printf("%lld\n", ans * inv[k + 1] % mod);}return 0;
}

考虑生成函数

$∑i=0nCn+1iBi=[n=0]对两边加上Bn+1∑i=0n+1Cn+1iBi=[n=0]+Bn+1∑i=0nCniBi=[n=1]+Bn∑i=0n1i!(n−i)!Bi=[n=1]+Bnn!设指数型生成函数B(x)=∑n≥0Bnxn有B(x)ex=x+B(x)从而有B(x)=xex−1B(x)=x∑n≥01n!xn−1x∑n≥11n!xn1∑n≥01(n+1)!xn\sum_{i = 0} ^{n} C_{n + 1} ^{i} B_i = [n = 0]\\ 对两边加上B_{n + 1}\\ \sum_{i = 0} ^{n + 1} C_{n + 1} ^{i} B_i = [n = 0] + B_{n + 1}\\ \sum_{i = 0} ^{n} C_{n} ^{i} B_i = [n = 1] + B_{n}\\ \sum_{i = 0} ^{n} \frac{1}{i!(n - i)!} B_{i} = [n = 1] + \frac{B_{n}}{n!}\\ 设指数型生成函数B(x) = \sum_{n \geq 0} B_n x ^ n\\ 有B(x) e ^x = x + B(x)\\ 从而有B(x) = \frac{x}{e ^ x - 1}\\ B(x) = \frac{x}{\sum\limits_{n \geq 0} \frac{1}{n!} x ^ n - 1}\\ \frac{x}{\sum\limits_{n \geq 1} \frac{1}{n!} x ^ n}\\ \frac{1}{\sum\limits_{n \geq 0} \frac{1}{(n + 1)!}x ^ n}\\$
只需多项式求逆即可解决。

#include <bits/stdc++.h>using namespace std;const int N = 1e6 + 10, mod = 998244353;int r[N], b[N], t[N];int quick_pow(int a, int n) {int ans = 1;while (n) {if (n & 1) {ans = 1ll * a * ans % mod;}a = 1ll * a * a % mod;n >>= 1;}return ans;
}void get_r(int lim) {for (int i = 0; i < lim; i++) {r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);}
}void NTT(int *f, int lim, int rev) {for (int i = 0; i < lim; i++) {if (i < r[i]) {swap(f[i], f[r[i]]);}}for (int mid = 1; mid < lim; mid <<= 1) {int wn = quick_pow(3, (mod - 1) / (mid << 1));for (int len = mid << 1, cur = 0; cur < lim; cur += len) {int w = 1;for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;}}}if (rev == -1) {int inv = quick_pow(lim, mod - 2);reverse(f + 1, f + lim);for (int i = 0; i < lim; i++) {f[i] = 1ll * f[i] * inv % mod;}}
}void polyinv(int *f, int *g, int n) {if (n == 1) {g[0] = quick_pow(f[0], mod - 2);return ;}polyinv(f, g, n + 1 >> 1);for (int i = 0; i < n; i++) {t[i] = f[i];}int lim = 1;while (lim < 2 * n) {lim <<= 1;}get_r(lim);NTT(t, lim, 1);NTT(g, lim, 1);for (int i = 0; i < lim; i++) {int cur = (2 - 1ll * g[i] * t[i] % mod + mod) % mod;g[i] = 1ll * g[i] * cur % mod;t[i] = 0;}NTT(g, lim, -1);for (int i = n; i < lim; i++) {g[i] = 0;}
}int fac[N], ifac[N], B[N], n;void init() {fac[0] = 1;for (int i = 1; i < N; i++) {fac[i] = 1ll * fac[i - 1] * i % mod;}ifac[N - 1] = quick_pow(fac[N - 1], mod - 2);for (int i = N - 2; i >= 0; i--) {ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;}
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);init();n = 100000;for (int i = 0; i <= n + 1; i++) {B[i] = ifac[i + 1];}polyinv(B, b, n + 2);for (int i = 0; i <= n + 1; i++) {B[i] = 1ll * b[i] * fac[i] % mod;}for (int i = 0; i <= n + 1; i++) {printf("%d%c", B[i], i == n + 1 ? '\n' : ' ');}return 0;
}

P3711 仓鼠的数学题

$\sum_{k = 0} ^{n} S_k(x) a_k\\ 原本定义的S_k(x) = \sum_{i = 0} ^{x} i ^ k\\ 根据伯努利数的定义S'_k(x) = \sum_{i = 0} ^{x - 1} i ^ k\\ 则我们求F(x) = \sum_{k = 0} ^{n} S'_k(x) a_k, 答案即为F(x + 1)\\ 考虑先求F(x) = \sum_{k = 0} ^{n} a_k \frac{1}{k + 1} \sum\limits_{i = 0} ^{k}C_{k + 1} ^{i}B_i x ^{k - i + 1}\\ \sum_{k = 0} ^{n} \frac{a_k}{k + 1} \sum_{i = 0} ^{k} \frac{(k + 1)!}{i!(k + 1 - i)!} B_i x ^{k - i + 1}\\ \sum_{k = 0} ^{n} a_k k! \sum_{i = 0} ^{k} \frac{B_i}{i!}\frac{x ^{k - i + 1}}{(k - i + 1)!}\\ \sum_{i = 1} ^{n + 1} \frac{x ^{i}}{i!} \sum_{k = i - 1} ^{n} \left(a_k k!\right) \left( \frac{B_{k + 1 - i}}{(k + 1 - i)!} \right)\\ 设f(n) = a_k k!, g(n) = \frac{B_{n}}{n!}, 另g(n) = g(n - i)即翻转一下\\ \sum_{i = 1} ^{n + 1} \frac{x ^ i}{i!} \sum_{k = i - 1} ^{n} f(k) g(n + i - k)\\ 容易发现后面是一个卷积，所以可以优化\\$
下面的 $n = n + 1$ ，由上易知 $F (x)$ 是一个 $n + 1$ 次多项式。
$\sum_{i = 0} ^{n} a_i x ^ i\\ F(x + 1) = \sum_{i = 0} ^{n} a_i (x + 1) ^ i\\ \sum_{i = 0} ^{n} a_i \sum_{j = 0} ^{i} C_{j} ^{i} x ^ j\\ \sum_{i = 0} ^{n} a_i \sum_{j = 0} ^{i} \frac{i!}{j!(i - j)!} x ^ j\\ \sum_{i = 0} ^{n} \frac{x ^ i}{i!} \sum_{j = i} ^{n} a_j j! \frac{1}{(j - i)!}\\ 另f(n) = a_n n!, g(n) = \frac{1}{n!},再翻转g(n)\\ \sum_{i = 0} ^{n} \frac{x ^{i}}{i!} \sum_{j = i} ^{n} f(j) g(n + i - j)\\ 后面也同样是一个卷积，所以可以优化\\$

#include <bits/stdc++.h>using namespace std;const int N = 1e6 + 10, mod = 998244353;int r[N], t[N];int quick_pow(int a, int n) {int ans = 1;while (n) {if (n & 1) {ans = 1ll * a * ans % mod;}a = 1ll * a * a % mod;n >>= 1;}return ans;
}void get_r(int lim) {for (int i = 0; i < lim; i++) {r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);}
}void NTT(int *f, int lim, int rev) {for (int i = 0; i < lim; i++) {if (i < r[i]) {swap(f[i], f[r[i]]);}}for (int mid = 1; mid < lim; mid <<= 1) {int wn = quick_pow(3, (mod - 1) / (mid << 1));for (int len = mid << 1, cur = 0; cur < lim; cur += len) {int w = 1;for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;}}}if (rev == -1) {int inv = quick_pow(lim, mod - 2);reverse(f + 1, f + lim);for (int i = 0; i < lim; i++) {f[i] = 1ll * f[i] * inv % mod;}}
}void polyinv(int *f, int *g, int n) {if (n == 1) {g[0] = quick_pow(f[0], mod - 2);return ;}polyinv(f, g, n + 1 >> 1);for (int i = 0; i < n; i++) {t[i] = f[i];}int lim = 1;while (lim < 2 * n) {lim <<= 1;}get_r(lim);NTT(t, lim, 1);NTT(g, lim, 1);for (int i = 0; i < lim; i++) {int cur = (2 - 1ll * g[i] * t[i] % mod + mod) % mod;g[i] = 1ll * g[i] * cur % mod;t[i] = 0;}NTT(g, lim, -1);for (int i = n; i < lim; i++) {g[i] = 0;}
}int fac[N], ifac[N], a[N], f[N], g[N], B[N], n;void init() {fac[0] = 1;for (int i = 1; i < N; i++) {fac[i] = 1ll * fac[i - 1] * i % mod;}ifac[N - 1] = quick_pow(fac[N - 1], mod - 2);for (int i = N - 2; i >= 0; i--) {ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;}
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);init();scanf("%d", &n);for (int i = 0; i <= n; i++) {scanf("%d", &a[i]);}for (int i = 0; i <= n + 1; i++) {B[i] = ifac[i + 1];}polyinv(B, g, n + 2);for (int i = 0; i <= n; i++) {f[i] = 1ll * a[i] * fac[i] % mod;}reverse(g, g + n + 2);int lim = 1;while (lim < 2 * n + 10) {lim <<= 1;}get_r(lim);NTT(f, lim, 1), NTT(g, lim, 1);for (int i = 0; i < lim; i++) {f[i] = 1ll * f[i] * g[i] % mod;}NTT(f, lim, -1);for (int i = 1; i <= n + 1; i++) {a[i] = 1ll * ifac[i] * f[n + i] % mod;}a[0] = 0;n++;for (int i = 0; i < lim; i++) {f[i] = g[i] = 0;}for (int i = 0; i <= n; i++) {f[i] = 1ll * a[i] * fac[i] % mod;g[i] = ifac[n - i];}NTT(f, lim, 1), NTT(g, lim, 1);for (int i = 0; i < lim; i++) {f[i] = 1ll * f[i] * g[i] % mod;}NTT(f, lim, -1);for (int i = 0; i <= n; i++) {printf("%lld%c", 1ll * ifac[i] * f[n + i] % mod, i == n ? '\n' : ' ');}return 0;
}