瞎寫的理性愉悅：正整數冪和與伯努利數

小粉兔原本對數學方面的知識（數論，計數）比較感興趣，直到他考了 CTS2019 和 APIO2019，
發現他的數學居然連沒怎么學的數據結構都不如，於是他自閉了。

考慮到從現在開始 CNOI 可能會變得計數和推式子化，小粉兔很可能要涼掉了。
現在補救可能還是來得及的（？），所以他推了一些式子嘗試救起瀕死的數學水平。

注意：以下內容極其低端，僅能體現 OIer 中數學的最低水平，請謹慎閱讀！

注意：以下內容極其低端，僅能體現 OIer 中數學的最低水平，請謹慎閱讀！

注意：以下內容極其低端，僅能體現 OIer 中數學的最低水平，請謹慎閱讀！

---

考慮計算這樣一個東西：\(\displaystyle\sum_{i=1}^{n}i^k\)。

即前 \(n\) 個正整數冪和，其中 \(0\le k\le 10^6\)，\(1\le n\le 10^{18}\) 且均為整數。

設 \(\displaystyle S_k(n)=\sum_{i=1}^{n}i^k\)，則顯然對於 \(k\in\mathbb{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\right]\right\}\)。
• \(\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}\right]\right\}\)。
• \(\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}\right]\right\}\)。
• \(\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}\right]\right\}\)。

---

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

觀察到 \(S_k(n)\) 中包含許多 \(i^k\)，可以聯想到其出現在某些有相似特征的生成函數中，
比如 \(\displaystyle\mathbf{OGF}\left\{\left[1,a,a^2,a^3,\ldots\right]\right\}=\frac{1}{1-ax}\) 或 \(\displaystyle\mathbf{EGF}\left\{\left[1,a,a^2,a^3,\ldots\right]\right\}=e^{ax}\)。

我們不妨選用 \(\mathbf{EGF}:e^{ax}\)。

考慮 \(\displaystyle\hat{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\right]\right\}\)，則有如下式子：

\[\begin{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} \]

最后的 \(\displaystyle\sum_{i=1}^{n}e^{ix}\) 可以這樣化簡：

\[\begin{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\\\therefore\hat{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{\displaystyle\sum_{i=0}^{\infty}\frac{n^{i+1}}{(i+1)!}x^i}{\displaystyle\sum_{i=0}^{\infty}\frac{(-1)^{i+1}}{(i+1)!}x^i}\)。

據此可以在 \(\mathcal{O}(m\log m)\) 的時間內，使用多項式求逆和卷積，求出 \(n\) 固定時 \(0\le k<m\) 的所有 \(S_k(n)\)。

---

我們之前提到了 \(S_k(x)\) 是關於 \(x\) 的 \(k+1\) 次多項式，假設固定了 \(k\) 能否求出此多項式呢？

當然我們無法對 \(1\le k\le 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\right]\right\}\)，注意到 \(\hat{G}(n,x)\) 的分子就包含此式，順着思路化簡：

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

則右邊：\(\displaystyle\frac{e^{nx}-1}{x}=\mathbf{EGF}\left\{\left[n,\frac{n^2}{2},\frac{n^3}{3},\ldots\right]\right\}\)。

設左邊的生成函數對應的級數為 \(B\)，即：

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

我們知道指數型生成函數的乘積對應了其對應級數的二項卷積，我們把式子寫成二項卷積的形式：

\[\begin{aligned}\hat{G}(n,x)&=\hat{B}(x)\cdot\mathbf{EGF}\left\{\left[n,\frac{n^2}{2},\frac{n^3}{3},\ldots\right]\right\}\\\therefore\forall k\ge 0,S_k(n)&=\sum_{j=0}^{k}\binom{k}{j}B_j\frac{n^{k-j+1}}{k-j+1}\\\therefore\sum_{i=1}^{n}i^k&=\frac{1}{k+1}\sum_{j=0}^{k}\binom{k+1}{j}B_jn^{k-j+1}\end{aligned} \]

則得到最終的式子：\(\displaystyle\sum_{i=1}^{n}i^k=\frac{1}{k+1}\sum_{j=0}^{k}\binom{k+1}{j}B_jn^{k-j+1}\)。

其中 \(B\) 就是所稱的“伯努利數”：\(\displaystyle\mathbf{EGF}\left\{B\right\}=\frac{xe^x}{e^x-1}\)。

其前 \(15\) 項為：

\[\begin{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{\displaystyle\sum_{i=0}^{\infty}\frac{1}{i!}x^i}{\displaystyle\sum_{i=0}^{\infty}\frac{1}{(i+1)!}x^i}\)。

這是適宜使用多項式求逆和卷積在 \(\mathcal{O}(n\log n)\) 的時間復雜度內求出前 \(n\) 個伯努利數的方法。

---

想求伯努利數，不想寫多項式，怎么辦？

沒關系，我們可以通過生成函數反推出遞推式，考慮如下式子：

\[\begin{aligned}\hat{B}(x)\cdot\frac{e^x-1}{x}&=e^x\\\mathbf{EGF}\left\{B\right\}\cdot\mathbf{EGF}\left\{\left[1,\frac{1}{2},\frac{1}{3},\ldots\right]\right\}&=\mathbf{EGF}\left\{\left[1,1,1,\ldots\right]\right\}\\\sum_{i=0}^{n}\binom{n}{i}B_i\frac{1}{n-i+1}&=1\\\sum_{i=0}^{n}\binom{n+1}{i}B_i&=n+1\end{aligned} \]

那么以 \(B_0=1\) 為邊界，按照遞推式 \(\displaystyle B_n=1-\frac{1}{n+1}\sum_{i=0}^{n-1}\binom{n+1}{i}B_i\)，
即可依次計算出 \(B_i\) 的值，假設需要處理前 \(m\) 個伯努利數，則復雜度為 \(\mathcal{O}(m^2)\)。

---

參考資料：
Bernoulli number, Wikipedia。
康復計划#3 簡單常用的幾種計算自然數冪和的方法, MoebiusMeow。
@_rqy 在 Universal OJ 用戶群中的解答。