瞎写的理性愉悦：正整数幂和与伯努利数

小粉兔原本对数学方面的知识（数论，计数）比较感兴趣，直到他考了 CTS2019 和 APIO2019，
发现他的数学居然连没怎么学的数据结构都不如，于是他自闭了。

考虑到从现在开始 CNOI 可能会变得计数和推式子化，小粉兔很可能要凉掉了。
现在补救可能还是来得及的（？），所以他推了一些式子尝试救起濒死的数学水平。

注意：以下内容极其低端，仅能体现 OIer 中数学的最低水平，请谨慎阅读！

考虑计算这样一个东西：(displaystylesum_{i=1}^{n}i^k)。

即前 (n) 个正整数幂和，其中 (0le kle 10^6)，(1le nle 10^{18}) 且均为整数。

设 (displaystyle S_k(n)=sum_{i=1}^{n}i^k)，则显然对于 (kinmathbb{N})，(S_k(x)) 是关于 (x) 的一个 (k+1) 次多项式。

举几个例子：

(displaystyle S_0(n)=overset{n}{overbrace{1+1+cdots+1}}=n=mathbf{OGF}left{left[0,1 ight] ight})。
(displaystyle S_1(n)=1+2+cdots+n=frac{n(n+1)}{2}=frac{1}{2}n+frac{1}{2}n^2=mathbf{OGF}left{left[0,frac{1}{2},frac{1}{2} ight] ight})。
(displaystyle S_2(n)=1^2+2^2+cdots+n^2=frac{n(n+1)(2n+1)}{6}=frac{1}{6}n+frac{1}{2}n^2+frac{1}{3}n^3=mathbf{OGF}left{left[0,frac{1}{6},frac{1}{2},frac{1}{3} ight] ight})。
(displaystyle S_3(n)=1^3+2^3+cdots+n^3=frac{n^2(n+1)^2}{4}=frac{1}{4}n^2+frac{1}{2}n^3+frac{1}{4}n^4=mathbf{OGF}left{left[0,0,frac{1}{4},frac{1}{2},frac{1}{4} ight] ight})。

要如何求出 (S_k(n)) 呢？

观察到 (S_k(n)) 中包含许多 (i^k)，可以联想到其出现在某些有相似特征的生成函数中，
比如 (displaystylemathbf{OGF}left{left[1,a,a^2,a^3,ldots ight] ight}=frac{1}{1-ax}) 或 (displaystylemathbf{EGF}left{left[1,a,a^2,a^3,ldots ight] ight}=e^{ax})。

我们不妨选用 (mathbf{EGF}:e^{ax})。

考虑 (displaystylehat{G}(n,x)=sum_{k=0}^{infty}frac{S_{k}(n)}{k!}x^k)，即 (mathbf{EGF}left{left[S_0(n),S_1(n),S_2(n),ldots ight] ight})，则有如下式子：

[egin{aligned}hat{G}(n,x)&=sum_{k=0}^{infty}frac{S_{k}(n)}{k!}x^k\&=sum_{k=0}^{infty}frac{sum_{i=1}^{n}i^k}{k!}x^k\&=sum_{i=1}^{n}sum_{k=0}^{infty}frac{i^k}{k!}x^k\&=sum_{i=1}^{n}e^{ix}end{aligned} ]

最后的 (displaystylesum_{i=1}^{n}e^{ix}) 可以这样化简：

[egin{aligned}(1-e^{-x})hat{G}(n,x)&=sum_{i=1}^{n}e^{ix}-sum_{i=1}^{n}e^{(i-1)x}\&=sum_{i=1}^{n} e^{ix}-sum_{i=0}^{n-1}e^{ix}\&=e^{nx}-1\ hereforehat{G}(n,x)&=frac{e^{nx}-1}{1-e^{-x}}end{aligned} ]

那么 (S_k(n)) 等于 (k![x^k]hat{G}(n,x))，也就等于 (displaystyle k![x^k]frac{displaystylesum_{i=0}^{infty}frac{n^{i+1}}{(i+1)!}x^i}{displaystylesum_{i=0}^{infty}frac{(-1)^{i+1}}{(i+1)!}x^i})。

据此可以在 (mathcal{O}(mlog m)) 的时间内，使用多项式求逆和卷积，求出 (n) 固定时 (0le k<m) 的所有 (S_k(n))。

我们之前提到了 (S_k(x)) 是关于 (x) 的 (k+1) 次多项式，假设固定了 (k) 能否求出此多项式呢？

当然我们无法对 (1le kle m) 的所有 (k) 求出所有多项式的系数了，因为有 (mathcal{O}(m^2)) 个系数。

仍然考察 (hat{G}(n,x))，我们希望的是将其写作一个与 (n) 无关的生成函数与一个与 (n^i)有关的生成函数的乘积。

考虑 (e^{nx}=mathbf{EGF}left{left[1,n,n^2,n^3,ldots ight] ight})，注意到 (hat{G}(n,x)) 的分子就包含此式，顺着思路化简：

[egin{aligned}hat{G}(n,x)&=frac{e^{nx}-1}{1-e^{-x}}\&=left(frac{xe^x}{e^x-1} ight)left(frac{e^{nx}-1}{x} ight)end{aligned} ]

则右边：(displaystylefrac{e^{nx}-1}{x}=mathbf{EGF}left{left[n,frac{n^2}{2},frac{n^3}{3},ldots ight] ight})。

设左边的生成函数对应的级数为 (B)，即：

[frac{xe^x}{e^x-1}=hat{B}(x)=mathbf{EGF}left{left[B_0,B_1,B_2,ldots ight] ight}=sum_{i=0}^{infty}frac{B_i}{i!}x^i ]

我们知道指数型生成函数的乘积对应了其对应级数的二项卷积，我们把式子写成二项卷积的形式：

[egin{aligned}hat{G}(n,x)&=hat{B}(x)cdotmathbf{EGF}left{left[n,frac{n^2}{2},frac{n^3}{3},ldots ight] ight}\ hereforeforall kge 0,S_k(n)&=sum_{j=0}^{k}inom{k}{j}B_jfrac{n^{k-j+1}}{k-j+1}\ hereforesum_{i=1}^{n}i^k&=frac{1}{k+1}sum_{j=0}^{k}inom{k+1}{j}B_jn^{k-j+1}end{aligned} ]

则得到最终的式子：(displaystylesum_{i=1}^{n}i^k=frac{1}{k+1}sum_{j=0}^{k}inom{k+1}{j}B_jn^{k-j+1})。

其中 (B) 就是所称的“伯努利数”：(displaystylemathbf{EGF}left{B ight}=frac{xe^x}{e^x-1})。

其前 (15) 项为：

[egin{aligned}B_0&=1&B_1&=frac{1}{2}&B_2&=frac{1}{6}&B_3&=0&B_4&=-frac{1}{30}\B_5&=0&B_6&=frac{1}{42}&B_7&=0&B_8&=-frac{1}{30}&B_9&=0\B_{10}&=frac{5}{66}&B_{11}&=0&B_{12}&=-frac{691}{2730}&B_{13}&=0&B_{14}&=frac{7}{6}end{aligned} ]

容易看出除了 (displaystyle B_1=frac{1}{2})，其它奇数位置的项均为 (0)。

考虑如何求出伯努利数。

从生成函数的角度考虑，有：(displaystyle B_n=n![x^n]frac{displaystylesum_{i=0}^{infty}frac{1}{i!}x^i}{displaystylesum_{i=0}^{infty}frac{1}{(i+1)!}x^i})。

这是适宜使用多项式求逆和卷积在 (mathcal{O}(nlog n)) 的时间复杂度内求出前 (n) 个伯努利数的方法。

想求伯努利数，不想写多项式，怎么办？

没关系，我们可以通过生成函数反推出递推式，考虑如下式子：

[egin{aligned}hat{B}(x)cdotfrac{e^x-1}{x}&=e^x\mathbf{EGF}left{B ight}cdotmathbf{EGF}left{left[1,frac{1}{2},frac{1}{3},ldots ight] ight}&=mathbf{EGF}left{left[1,1,1,ldots ight] ight}\sum_{i=0}^{n}inom{n}{i}B_ifrac{1}{n-i+1}&=1\sum_{i=0}^{n}inom{n+1}{i}B_i&=n+1end{aligned} ]

那么以 (B_0=1) 为边界，按照递推式 (displaystyle B_n=1-frac{1}{n+1}sum_{i=0}^{n-1}inom{n+1}{i}B_i)，
即可依次计算出 (B_i) 的值，假设需要处理前 (m) 个伯努利数，则复杂度为 (mathcal{O}(m^2))。

参考资料：
Bernoulli number, Wikipedia。
康复计划#3 简单常用的几种计算自然数幂和的方法, MoebiusMeow。
@_rqy 在 Universal OJ 用户群中的解答。