O(m²) 实现等幂求和——伯努利数

等幂求和

伯努利数是以雅各布

⋅

\cdot

⋅ 伯努利的名字命名的，我们记

S

m

(

n

)

=

∑

k

=

0

n

−

1

k

m

S_m(n) = \sum_{k = 0}^{n – 1} k^m

Sm​(n)=k=0∑n−1​km
伯努利 暴力 算出了前几项，并进行了观察：

S

0

(

n

)

=

n

S

1

(

n

)

=

1

2

n

2

−

1

2

n

S

2

(

n

)

=

1

3

n

3

−

1

2

n

2

+

1

6

n

S

3

(

n

)

=

1

4

n

4

−

1

2

n

3

+

1

4

n

2

+

0

n

S

4

(

n

)

=

1

5

n

5

−

1

2

n

4

+

1

3

n

3

+

0

n

2

−

1

30

n

S

5

(

n

)

=

1

6

n

6

−

1

2

n

5

+

5

12

n

4

+

0

n

3

−

1

12

n

2

+

0

n

S

6

(

n

)

=

1

7

n

7

−

1

2

n

6

+

1

2

n

5

+

0

n

4

−

1

6

n

3

+

0

n

2

+

1

42

n

\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}

S0​(n)=nS1​(n)=21​n2−21​nS2​(n)=31​n3−21​n2+61​nS3​(n)=41​n4−21​n3+41​n2+0nS4​(n)=51​n5−21​n4+31​n3+0n2−301​nS5​(n)=61​n6−21​n5+125​n4+0n3−121​n2+0nS6​(n)=71​n7−21​n6+21​n5+0n4−61​n3+0n2+421​n​
我们发现，在

S

m

(

n

)

S_m(n)

Sm​(n) 中，

n

m

+

1

n^{m + 1}

nm+1 的系数总是

1

m

+

1

\dfrac{1}{m +1}

m+11​，

n

m

n^m

nm 的系数是

−

1

2

-\dfrac{1}{2}

−21​，

n

m

−

1

n^{m – 1}

nm−1 的系数是

m

12

\dfrac{m}{12}

12m​，

n

m

−

2

n^{m – 2}

nm−2 的系数是

0

0

0，

n

m

−

3

n^{m – 3}

nm−3 的系数是

−

m

(

m

−

1

)

(

m

−

2

)

720

– \dfrac{m (m – 1) (m – 2)}{720}

−720m(m−1)(m−2)​，

n

m

−

4

n^{m – 4}

nm−4 的系数是

0

⋯

⋯

0\cdots\cdots

0⋯⋯ 而

n

m

−

k

n^{m – k}

nm−k 的系数总是某个常数乘上

m

k

‾

=

m

!

(

m

−

k

)

!

m^{\underline{k}} = \dfrac{m!}{(m – k)!}

mk​=(m−k)!m!​。

递推公式

根据伯努利的发现（值得一提的是，他并没有给出证明），我们有

S

m

(

n

)

=

1

m

+

1

∑

k

=

0

m

(

m

+

1

k

)

B

k

n

m

+

1

−

k

S_m(n) = \dfrac{1}{m + 1} \sum_{k = 0}^m \dbinom{m + 1}{k} B_k n^{m + 1 – k}

Sm​(n)=m+11​k=0∑m​(km+1​)Bk​nm+1−k
其中

B

n

B_n

Bn​ 为伯努利数，其由一个递归关系定义：

∑

k

=

0

m

(

m

+

1

k

)

B

k

=

[

m

=

0

]

\sum_{k = 0}^m \dbinom{m + 1}{k} B_k = [m = 0]

k=0∑m​(km+1​)Bk​=[m=0]
其实你也可以写成

B

0

=

1

∑

k

=

0

m

(

m

+

1

k

)

B

k

=

0

,

m

≥

1

B_0 = 1 \\ \sum_{k = 0}^m \dbinom{m + 1}{k} B_k = 0, m \ge 1

B0​=1k=0∑m​(km+1​)Bk​=0,m≥1
这样方便你进行计算。

根据手算，不难得出前几个值：

$n$	$0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$	$10$	$11$	$12$
$B_n$	$1$	$-\frac{1}{2}$	$\frac{1}{6}$	$0$	$-\frac{1}{30}$	$0$	$\frac{1}{42}$	$0$	$-\frac{1}{30}$	$0$	$\frac{5}{66}$	$0$	$-\frac{691}{2730}$

（其中

B

12

B_{12}

B12​ 这个奇怪分数的出现真是给那些有关

B

n

B_n

Bn​ 的简单封闭形式的猜想泼冷水。）

证明：

首先令

S

^

m

(

n

)

=

1

m

+

1

∑

k

=

0

m

(

m

+

1

k

)

B

k

n

m

+

1

−

k

\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)=m+11​k=0∑m​(km+1​)Bk​nm+1−k
那么现在我们就是要证明

S

m

(

n

)

=

S

^

m

(

n

)

S_m(n) = \hat{S}_m(n)

Sm​(n)=S^m​(n)，考虑使用数学归纳法。

首先

S

0

(

n

)

=

∑

i

=

0

n

−

1

n

0

=

n

S

^

0

(

n

)

=

1

1

(

1

0

)

B

0

n

1

=

n

\begin{array}{ll} S_0(n) = \sum\limits_{i = 0}^{n – 1} n^0 = n \\ \hat{S}_0(n) = \dfrac{1}{1} \dbinom{1}{0} B_0 n^{1} = n \end{array}

S0​(n)=i=0∑n−1​n0=nS^0​(n)=11​(01​)B0​n1=n​

S

m

+

1

(

n

)

+

n

m

+

1

=

∑

k

=

0

n

k

m

+

1

=

∑

k

=

0

n

−

1

(

k

+

1

)

m

+

1

=

∑

k

=

0

n

−

1

∑

j

=

0

m

+

1

(

m

+

1

j

)

k

j

=

∑

j

=

0

m

+

1

(

m

+

1

j

)

∑

k

=

0

n

−

1

k

j

=

∑

j

=

0

m

+

1

(

m

+

1

j

)

S

j

(

n

)

\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}

Sm+1​(n)+nm+1​=k=0∑n​km+1=k=0∑n−1​(k+1)m+1=k=0∑n−1​j=0∑m+1​(jm+1​)kj=j=0∑m+1​(jm+1​)k=0∑n−1​kj=j=0∑m+1​(jm+1​)Sj​(n)​

当

m

≥

1

m\ge 1

m≥1 时，假设

∀

i

∈

[

0

,

m

)

\forall i \in [0, m)

∀i∈[0,m) 均有

S

i

(

n

)

=

S

^

i

(

n

)

S_i(n) = \hat{S}_i(n)

Si​(n)=S^i​(n)，将

S

m

+

1

(

n

)

S_{m +1}(n)

Sm+1​(n) 移项：

n

m

+

1

=

∑

j

=

0

m

+

1

(

m

+

1

j

)

S

j

(

n

)

−

S

m

+

1

(

n

)

=

∑

j

=

0

m

(

m

+

1

j

)

S

j

(

n

)

=

∑

j

=

0

m

−

1

(

m

+

1

j

)

S

j

(

n

)

+

(

m

+

1

m

)

S

m

(

n

)

=

∑

j

=

0

m

−

1

(

m

+

1

j

)

S

^

j

(

n

)

+

(

m

+

1

)

S

m

(

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}

nm+1​=j=0∑m+1​(jm+1​)Sj​(n)−Sm+1​(n)=j=0∑m​(jm+1​)Sj​(n)=j=0∑m−1​(jm+1​)Sj​(n)+(mm+1​)Sm​(n)=j=0∑m−1​(jm+1​)S^j​(n)+(m+1)Sm​(n)​
设

Δ

=

S

m

(

n

)

−

S

^

m

(

n

)

\Delta = S_m(n) – \hat{S}_m(n)

Δ=Sm​(n)−S^m​(n)
那么现在就是要证

Δ

=

0

\Delta = 0

Δ=0。

n

m

+

1

=

[

∑

j

=

0

m

−

1

(

m

+

1

j

)

S

^

j

(

n

)

+

(

m

+

1

)

S

^

m

(

n

)

]

+

(

m

+

1

)

S

m

(

n

)

−

(

m

+

1

)

S

^

m

(

n

)

=

∑

j

=

0

m

(

m

+

1

j

)

S

^

j

(

n

)

+

(

m

+

1

)

Δ

\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}

nm+1​=[j=0∑m−1​(jm+1​)S^j​(n)+(m+1)S^m​(n)]+(m+1)Sm​(n)−(m+1)S^m​(n)=j=0∑m​(jm+1​)S^j​(n)+(m+1)Δ​
根据

S

^

m

(

n

)

\hat{S}_m(n)

S^m​(n) 的定义将其展开，有

n

m

+

1

=

∑

j

=

0

m

(

m

+

1

j

)

1

j

+

1

∑

k

=

0

j

(

j

+

1

k

)

B

k

n

j

+

1

−

k

+

(

m

+

1

)

Δ

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

nm+1=j=0∑m​(jm+1​)j+11​k=0∑j​(kj+1​)Bk​nj+1−k+(m+1)Δ
将

k

k

k 改为倒序枚举 并 利用 对称恒等式

(

n

m

)

=

(

n

n

−

m

)

,

n

∈

N

,

k

∈

Z

\dbinom{n}{m} = \dbinom{n}{n – m}, n\in \mathbb{N}, k \in \mathbb{Z}

(mn​)=(n−mn​),n∈N,k∈Z

n

m

+

1

=

∑

j

=

0

m

(

m

+

1

j

)

1

j

+

1

∑

k

=

0

j

(

j

+

1

j

−

k

)

B

j

−

k

n

k

+

1

+

(

m

+

1

)

Δ

=

∑

j

=

0

m

(

m

+

1

j

)

1

j

+

1

∑

k

=

0

j

(

j

+

1

k

+

1

)

B

j

−

k

n

k

+

1

+

(

m

+

1

)

Δ

\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}

nm+1​=j=0∑m​(jm+1​)j+11​k=0∑j​(j−kj+1​)Bj−k​nk+1+(m+1)Δ=j=0∑m​(jm+1​)j+11​k=0∑j​(k+1j+1​)Bj−k​nk+1+(m+1)Δ​
利用 吸收恒等式

(

r

k

)

=

r

k

(

r

−

1

k

−

1

)

,

k

∈

Z

,

k

≠

0

\dbinom{r}{k} = \dfrac{r}{k} \dbinom{r – 1}{k – 1}, k \in \mathbb{Z}, k \ne 0

(kr​)=kr​(k−1r−1​),k∈Z,k​=0

n

m

+

1

=

∑

j

=

0

m

(

m

+

1

j

)

1

j

+

1

∑

k

=

0

j

j

+

1

k

+

1

(

j

k

)

B

j

−

k

n

k

+

1

+

(

m

+

1

)

Δ

=

∑

j

=

0

m

(

m

+

1

j

)

∑

k

=

0

j

n

k

+

1

k

+

1

(

j

k

)

B

j

−

k

+

(

m

+

1

)

Δ

=

∑

k

=

0

m

n

k

+

1

k

+

1

∑

j

=

k

m

(

m

+

1

j

)

(

j

k

)

B

j

−

k

+

(

m

+

1

)

Δ

\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}

nm+1​=j=0∑m​(jm+1​)j+11​k=0∑j​k+1j+1​(kj​)Bj−k​nk+1+(m+1)Δ=j=0∑m​(jm+1​)k=0∑j​k+1nk+1​(kj​)Bj−k​+(m+1)Δ=k=0∑m​k+1nk+1​j=k∑m​(jm+1​)(kj​)Bj−k​+(m+1)Δ​

利用 三项式版恒等式

(

r

m

)

(

m

k

)

=

(

r

k

)

(

r

−

k

m

−

k

)

\dbinom{r}{m} \dbinom{m}{k} = \dbinom{r}{k} \dbinom{r – k}{m – k}

(mr​)(km​)=(kr​)(m−kr−k​)

n

m

+

1

=

∑

k

=

0

m

n

k

+

1

k

+

1

∑

j

=

k

m

(

m

+

1

k

)

(

m

−

k

+

1

j

−

k

)

B

j

−

k

+

(

m

+

1

)

Δ

=

∑

k

=

0

m

n

k

+

1

k

+

1

(

m

+

1

k

)

∑

j

=

k

m

(

m

−

k

+

1

j

−

k

)

B

j

−

k

+

(

m

+

1

)

Δ

\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}

nm+1​=k=0∑m​k+1nk+1​j=k∑m​(km+1​)(j−km−k+1​)Bj−k​+(m+1)Δ=k=0∑m​k+1nk+1​(km+1​)j=k∑m​(j−km−k+1​)Bj−k​+(m+1)Δ​

将

j

−

k

j – k

j−k 改为

j

j

j

n

m

+

1

=

∑

k

=

0

m

n

k

+

1

k

+

1

(

m

+

1

k

)

∑

j

=

0

m

−

k

(

m

−

k

+

1

j

)

B

j

+

(

m

+

1

)

Δ

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

nm+1=k=0∑m​k+1nk+1​(km+1​)j=0∑m−k​(jm−k+1​)Bj​+(m+1)Δ
又根据伯努利数的递归定义

∑

k

=

0

m

(

m

+

1

k

)

B

k

=

[

m

=

0

]

\sum_{k = 0}^m \dbinom{m + 1}{k} B_k = [m = 0]

k=0∑m​(km+1​)Bk​=[m=0]

n

m

+

1

=

∑

k

=

0

m

n

k

+

1

k

+

1

(

m

+

1

k

)

[

m

−

k

=

0

]

+

(

m

+

1

)

Δ

=

n

m

+

1

m

+

1

(

m

+

1

m

)

+

(

m

+

1

)

Δ

=

n

m

+

1

+

(

m

+

1

)

Δ

\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}

nm+1​=k=0∑m​k+1nk+1​(km+1​)[m−k=0]+(m+1)Δ=m+1nm+1​(mm+1​)+(m+1)Δ=nm+1+(m+1)Δ​

啊哈！至此，我们推出了

(

m

+

1

)

Δ

=

0

∵

m

+

1

≠

0

∴

Δ

=

0

(m + 1) \Delta = 0 \\ \because m + 1 \ne 0 \\ \therefore \Delta = 0

(m+1)Δ=0∵m+1​=0∴Δ=0
证完喽！

当然，还有一种用 指数生成函数 的更简单的证明方法，因为笔者才疏学浅不会生成函数就先不证了，到时候在生成函数的笔记里证（

代码实现

显然这东西一般用在有模数的情况下。

下列代码计算的是

∑

k

=

0

n

k

m

\sum\limits_{k = 0}^n k^m

k=0∑n​km，因此需要 n++;。

时间复杂度为

O

(

m

2

)

\Omicron(m^2)

O(m2)。

暴力计算的时间为

O

(

n

log

⁡

m

)

\Omicron(n\log m)

O(nlogm)，在

n

n

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; i++)
	{
		C[i][0] = C[i][i] = 1;
	}
	for (int i = 1; i <= N; 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.