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.