O(m²) 实现等幂求和——伯努利数
等幂求和
伯努利数是以雅各布 \(\cdot\) 伯努利的名字命名的，我们记
\[S_m(n) = \sum_{k = 0}^{n - 1} k^m
\]
伯努利 暴力 算出了前几项，并进行了观察：
\[\begin{array}{ll}
S_0(n) = n \\
S_1(n) = \dfrac{1}{2} n^2 - \dfrac{1}{2} n \\
S_2(n) = \dfrac{1}{3} n^3 - \dfrac{1}{2} n^2 + \dfrac{1}{6} n \\
S_3(n) = \dfrac{1}{4} n^4 - \dfrac{1}{2} n^3 + \dfrac{1}{4} n^2 + 0 n \\
S_4(n) = \dfrac{1}{5} n^5 - \dfrac{1}{2} n^4 + \dfrac{1}{3} n^3 + 0 n^2 - \dfrac{1}{30} n \\
S_5(n) = \dfrac{1}{6} n^6 - \dfrac{1}{2} n^5 + \dfrac{5}{12} n^4 + 0 n^3 - \dfrac{1}{12} n^2 + 0n \\
S_6(n) = \dfrac{1}{7} n^7 - \dfrac{1}{2} n^6 + \dfrac{1}{2} n^5 + 0 n^4 - \dfrac{1}{6} n^3 + 0 n^2 + \dfrac{1}{42} n
\end{array}
\]
我们发现，在 \(S_m(n)\) 中，\(n^{m + 1}\) 的系数总是 \(\dfrac{1}{m +1}\)，\(n^m\) 的系数是 \(-\dfrac{1}{2}\)，\(n^{m - 1}\) 的系数是 \(\dfrac{m}{12}\)，\(n^{m - 2}\) 的系数是 \(0\)，\(n^{m - 3}\) 的系数是 \(- \dfrac{m (m - 1) (m - 2)}{720}\)，\(n^{m - 4}\) 的系数是 \(0\cdots\cdots\) 而 \(n^{m - k}\) 的系数总是某个常数乘上 \(m^{\underline{k}} = \dfrac{m!}{(m - k)!}\)。
递推公式
根据伯努利的发现，我们有
\[S_m(n) = \dfrac{1}{m + 1} \sum_{k = 0}^m \dbinom{m + 1}{k} B_k n^{m + 1 - k}
\]
其中 \(B_n\) 为伯努利数，其由一个递归关系定义：
\[\sum_{k = 0}^m \dbinom{m + 1}{k} B_k = [m = 0]
\]
便于计算的表示：
\[B_0 = 1 \\
B_m = \frac{-\sum\limits_{k = 0}^{m - 1} \dbinom{m + 1}{k} B_k}{m + 1}, m > 0
\]
根据手算，不难得出前几个值：
| \(n\) | \(0\) | \(1\) | \(2\) | \(3\) | \(4\) | \(5\) | \(6\) | \(7\) | \(8\) | \(9\) | \(10\) | \(11\) | \(12\) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| \(B_n\) | \(1\) | \(-\dfrac{1}{2}\) | \(\dfrac{1}{6}\) | \(0\) | \(-\dfrac{1}{30}\) | \(0\) | \(\dfrac{1}{42}\) | \(0\) | \(-\dfrac{1}{30}\) | \(0\) | \(\dfrac{5}{66}\) | \(0\) | \(-\dfrac{691}{2730}\) |
证明：
首先令
\[\hat{S}_m(n) = \dfrac{1}{m + 1} \sum_{k = 0}^m \dbinom{m + 1}{k} B_k n^{m + 1 - k} \]那么现在我们就是要证明 \(S_m(n) = \hat{S}_m(n)\)，考虑使用数学归纳法。
首先
\[\begin{array}{ll} S_0(n) = \sum\limits_{i = 0}^{n - 1} i^0 = n \\ \hat{S}_0(n) = \dfrac{1}{1} \dbinom{1}{0} B_0 n^{1} = n \end{array} \]\[\begin{aligned} S_{m + 1}(n) + n^{m + 1} & = \sum_{k = 0}^n k^{m + 1} \\ & = \sum_{k = 0}^{n - 1} (k + 1)^{m + 1} \\ & = \sum_{k = 0}^{n - 1} \sum_{j = 0}^{m + 1} \dbinom{m + 1}{j} k^j \\ & = \sum_{j = 0}^{m + 1} \dbinom{m + 1}{j} \sum_{k = 0}^{n - 1} k^j \\ & = \sum_{j = 0}^{m + 1} \dbinom{m + 1}{j} S_j(n) \end{aligned} \]当 \(m\ge 1\) 时，假设 \(\forall i \in [0, m)\) 均有 \(S_i(n) = \hat{S}_i(n)\)，将 \(S_{m +1}(n)\) 移项：
\[\begin{aligned} n^{m + 1} & = \sum_{j = 0}^{m + 1} \dbinom{m + 1}{j} S_j(n) - S_{m + 1}(n) \\ & = \sum_{j = 0}^m \dbinom{m + 1}{j} S_j(n) \\ & = \sum_{j = 0}^{m - 1} \dbinom{m + 1}{j} S_j(n) + \dbinom{m + 1}{m} S_m(n) \\ & = \sum_{j = 0}^{m - 1} \dbinom{m + 1}{j} \hat{S}_j(n) + (m + 1) S_m(n) \end{aligned} \]设
\[\Delta = S_m(n) - \hat{S}_m(n) \]那么现在就是要证 \(\Delta = 0\)。
\[\begin{aligned} n^{m + 1} & = \left[\sum_{j = 0}^{m - 1} \dbinom{m + 1}{j} \hat{S}_j(n) + (m + 1) \hat{S}_m(n) \right] + (m + 1) S_m(n) - (m + 1) \hat{S}_m(n) \\ & = \sum_{j = 0}^m \dbinom{m + 1}{j} \hat{S}_j(n) + (m + 1) \Delta \end{aligned} \]根据 \(\hat{S}_m(n)\) 的定义将其展开，有
\[n^{m + 1} = \sum_{j = 0}^m \dbinom{m +1}{j} \dfrac{1}{j + 1} \sum_{k = 0}^j \dbinom{j + 1}{k} B_k n^{j + 1 - k} + (m + 1) \Delta \]将 \(k\) 改为倒序枚举 并 利用 对称恒等式
\[\dbinom{n}{m} = \dbinom{n}{n - m}, n\in \mathbb{N}, k \in \mathbb{Z} \]\[\begin{aligned} n^{m +1} & = \sum_{j = 0}^m \dbinom{m + 1}{j} \dfrac{1}{j + 1} \sum_{k = 0}^j \dbinom{j + 1}{j - k} B_{j - k} n^{k + 1} + (m + 1) \Delta \\ & = \sum_{j = 0}^m \dbinom{m + 1}{j} \dfrac{1}{j + 1} \sum_{k = 0}^j \dbinom{j + 1}{k + 1} B_{j - k} n^{k + 1} + (m + 1) \Delta \end{aligned} \]利用 吸收恒等式
\[\dbinom{r}{k} = \dfrac{r}{k} \dbinom{r - 1}{k - 1}, k \in \mathbb{Z}, k \ne 0 \]\[\begin{aligned} n^{m + 1} & = \sum_{j = 0}^m \dbinom{m + 1}{j} \dfrac{1}{j + 1} \sum_{k = 0}^j \dfrac{j + 1}{k + 1} \dbinom{j}{k} B_{j - k} n^{k + 1} + (m +1) \Delta \\ & = \sum_{j = 0}^m \dbinom{m + 1}{j} \sum_{k = 0}^j \dfrac{n^{k + 1}}{k + 1} \dbinom{j}{k} B_{j - k} + (m + 1) \Delta \\ & = \sum_{k = 0}^m \dfrac{n^{k + 1}}{k + 1} \sum_{j = k}^m \dbinom{m + 1}{j} \dbinom{j}{k} B_{j - k} + (m + 1) \Delta \end{aligned} \]利用 三项式版恒等式
\[\dbinom{r}{m} \dbinom{m}{k} = \dbinom{r}{k} \dbinom{r - k}{m - k} \]\[\begin{aligned} n^{m + 1} & = \sum_{k = 0}^m \dfrac{n^{k + 1}}{k + 1} \sum_{j = k}^m \dbinom{m + 1}{k} \dbinom{m - k + 1}{j - k} B_{j - k} + (m + 1) \Delta \\ & = \sum_{k = 0}^m \dfrac{n^{k + 1}}{k + 1} \dbinom{m + 1}{k} \sum_{j = k}^m \dbinom{m - k + 1}{j - k} B_{j - k} + (m + 1) \Delta \end{aligned} \]将 \(j - k\) 改为 \(j\)
\[n^{m + 1} = \sum_{k = 0}^m \dfrac{n^{k + 1}}{k + 1} \dbinom{m + 1}{k} \sum_{j = 0}^{m - k} \dbinom{m - k + 1}{j} B_j + (m + 1) \Delta \]又根据伯努利数的递归定义
\[\sum_{k = 0}^m \dbinom{m + 1}{k} B_k = [m = 0] \]\[\begin{aligned} n^{m + 1} & = \sum_{k = 0}^m \dfrac{n^{k + 1}}{k + 1} \dbinom{m + 1}{k} [m - k = 0] + (m + 1) \Delta \\ & = \dfrac{n^{m + 1}}{m + 1} \dbinom{m + 1}{m} + (m + 1) \Delta \\ & = n^{m + 1} + (m + 1) \Delta \end{aligned} \]至此，我们推出了
\[(m + 1) \Delta = 0 \\ \because m + 1 \ne 0 \\ \therefore \Delta = 0 \]证毕。
当然，还有一种用 指数生成函数 的更简单的证明方法，就先不证了（
代码实现
显然这东西一般用在有模数的情况下。
下列代码计算的是 \(\sum\limits_{k = 0}^n k^m\)，因此需要 n++;。
时间复杂度为 \(\Omicron(m^2)\)。
暴力计算的时间为 \(\Omicron(n\log m)\)，在 \(n\) 巨大时伯努利数有巨大优势，参见 P6271 [湖北省队互测2014]一个人的数论 使用莫反 + 伯努利数或时间复杂度同级的拉格朗日插值。
//18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;
const int MAXN = 1e4 + 5;
const int N = 1e4;
const int MOD = 1e9 + 7;
int add(int a, int b) {return (a + b) % MOD;}
int mul(int a, int b) {return (ll)a * b % MOD;}
int inv[MAXN], C[MAXN][MAXN], B[MAXN];
void init()
{
for (int i = 0; i <= N + 1; i++)
{
C[i][0] = C[i][i] = 1;
}
for (int i = 1; i <= N + 1; i++)
{
for (int j = 1; j < i; j++)
{
C[i][j] = add(C[i - 1][j], C[i - 1][j - 1]);
}
}
inv[1] = 1;
for (int i = 2; i <= N + 1; i++)
{
inv[i] = mul(MOD - MOD / i, inv[MOD % i]);
}
B[0] = 1;
for (int m = 1; m <= N; m++)
{
int sum = 0;
for (int k = 0; k < m; k++)
{
sum = add(sum, mul(C[m + 1][k], B[k]));
}
B[m] = mul(MOD - sum, inv[m + 1]);
}
}
int qpow(int a, int b)
{
int base = a, ans = 1;
while (b)
{
if (b & 1)
{
ans = mul(ans, base);
}
base = mul(base, base);
b >>= 1;
}
return ans;
}
int main()
{
init(); // 计算伯努利数
int n, m;
scanf("%d%d", &n, &m);
n++;
int res = 0;
for (int k = 0; k <= m; k++)
{
res += mul(mul(C[m + 1][k], B[k]), qpow(n, m + 1 - k));
}
printf("%d\n", mul(inv[m + 1], res));
return 0;
}
参考资料
- [1] 伯努利数. OI-wiki
- [2] Concrete Mathematics 5.1 BASIC IDENTITIES, 6.5 BERNOULLI NUMBER. Ronald L. Graham, Donald E. Knuth, Oren Patashnik.
