本篇口胡寫給我自己這樣的東西都忘光的殘廢選手 以及暫時還不會自然數冪和的人…
這里大概給出最簡單的幾種方法：擾動法(化為遞推式)，斯特林數(離散微積分)，高階差分(牛頓級數)，伯努利數(指數生成函數)…
不同方法的思維難度、普適程度、實現難度、時間復雜度上面都有差異…同時自然數冪和是探究各種求和方法的經典例子，了解多一點它的做法對於處理各種求和問題是有所幫助的…
問題：求$\sum_{k=0}^{n} k^t$，其中$t \in \mathbb{N}$是一個常數。要求求解的時間復雜度與$n$無關。（特別地，$0^0=1$）
方法0、 暴力
我們for一遍$0$到$n$，然后for一遍$0$到$t$就能算出來啦！比較慢？我們可以用快速冪算$k^t$，就不用for一次啦！時間復雜度$O\left( n \log t \right)$。
以上純屬賣萌。
這個復雜度是不能接受的，因為它與$n$有關…我們的要求就是與$n$無關，因為實際問題中$n$可以很大（弄幾百位的大數都可以），如果直接帶上一個$n$的話是無法快速求解的。
方法1、擾動法
名字可能有點不是那么眼熟，但是其實應該是老朋友了…也可以叫“寫兩遍”或者別的什么…處理等比數列的時候就和錯位相減類似（不過思維方向不同，所以擾動法更加普適）。
大概的思路是，把求和式寫兩遍放在左右兩邊，其中一邊提出末項，另一邊提出首項，然后嘗試用一邊表示另一邊（但是並不是總是能成功）。
令$S(n)=\sum_{k=0}^{n}a_k$，得到$S(n)+a_{n+1}=a_0+\sum_{k=1}^{n+1}a_k$。
這個方法可以直接處理等比數列。現在我們可以嘗試對自然數冪和使用這種辦法：
令$S_t(n)=\sum_{k=0}^{n}k^t$，得到$S_t(n)+(n+1)^t=0+\sum_{k=1}^{n+1}k^t$。
我們嘗試用$S_t(n)$表示右邊，所以試圖把右邊化為類似的形式。
$$ S_t(n)+(n+1)^t = \sum_{k=0}^{n}(k+1)^t $$
對於這種形式顯然可以使用二項式定理，右邊轉化成$\sum_{k=0}^{n} \sum_{i=0}^{t} \binom{t}{i} k^i $。
已經很類似了，接着我們交互求和次序：
$$\begin{aligned} S_t(n)+(n+1)^t & = \sum_{i=0}^{t} \binom{t}{i} \sum_{k=0}^{n} k^i \\ & = \sum_{i=0}^{t} \binom{t}{i} S_i(n) \end{aligned} $$
雖然和我們一開始想的有點出入，我們想用$S_t(n)$表示右邊，結果右邊的$S_t(n)$系數為$1$，所以它們消掉了…但是過程中我們引入了$S_i(n), i \leq t$，這意味着我們實際上找到了一個隱含的遞推式，我們可以用自然數的$0$到$i-2$次方和推出$i-1$次方和。（如果我們一開始嘗試求$S_{t+1}(n)$的話，我們得到的就是$S_t$的式子，不過這都是一樣，因為我們已經得到遞推式了，只是重新命名變量的區別）。所以我們用$t+1$代替式子中的$t$就可以得到：
$$ (n+1)^{t+1} = \binom{t+1}{t} S_t(n) + \sum_{i=0}^{t-1} \binom{t+1}{i} S_i(n) $$
移項得到：
$$ S_t(n) = \frac{1}{t+1} \left( (n+1)^{t+1} - \sum_{i=0}^{t-1} \binom{t+1}{i} S_i(n) \right) $$
這便是一個我們可以利用的遞推式了。我們知道$S_0(n)$就是$0$到$n$的自然數個數$n+1$，我們甚至可以不用知道$S_1(n)$的結論（雖然我們已經知道那是$n(n+1)/2$），因為我們可以遞推得到。所以對於一個確定的$n$，我們從$S_0(n)$開始一直遞推到$S_t(n)$就可以了。這樣子我們需要的時間是$O\left(t^2\right)$的，因為每項遞推都需要用到前面所有項的信息。這其實挺不錯了，因為它與$n$無關，而且也不是太大。
在后面的推導中我們可以發現自然數冪和其實是一個關於$n$的$t+1$次多項式。我們這個做法對於求多項式的系數其實不是很優的，需要$O\left(t^3\right)$的時間（因為不再是代一個數算而是代一個多項式算了）。如果知道這是一個多項式的話暴力高斯消元求系數甚至都可以做到$O\left( t^3 \right)$。所以我們來看下面這個做法。
方法2、斯特林數(離散微積分)
你問那個括號里是什么意思？其實就是對於方法的補充說明啦…
雖然一開始看起來比較高端，但是這個方法是所有方法里面最無腦的了…非常簡單，只要知道結論就會了…
首先簡單說一下斯特林數…斯特林數有兩類，第一類是輪換數，代表$n$個數排成$k$個非空的環的方案數，我們用$\left[ \begin{matrix}n\\k\end{matrix} \right]$表示；第二類是子集數，代表$n$個數形成$k$個非空集合的方案數，我們用$\left\{ \begin{matrix}n\\k\end{matrix} \right\}$表示。
根據意義我們可以得到它們的遞推式和邊界條件。顯然如果$n<k$，它們都是$0$；如果$n=k$，它們都是$1$；如果$n>0$且$k=0$，它們都是$0$。
加入一個新數的時候，可以考慮自己形成一個新的環，或者插入任意一個已有的數后面，所以
$$\left[ \begin{matrix}n\\k\end{matrix} \right]=\left[ \begin{matrix}n-1\\k-1\end{matrix} \right]+(n-1)\left[ \begin{matrix}n-1\\k\end{matrix} \right]$$
同樣，加入一個新數的時候，可以考慮自己形成一個新的集合，或者插入任意一個已有的集合。
$$\left\{ \begin{matrix}n\\k\end{matrix} \right\}=\left\{ \begin{matrix}n-1\\k-1\end{matrix} \right\}+k\left\{ \begin{matrix}n-1\\k\end{matrix} \right\}$$
我這次要用它們在我們使用的通常冪$n^k$和更加易於求和的下降冪$n^{\underline{k}}$之間轉化。其中$n^{\underline{k}}=n(n-1)(n-2) \cdots (n-k+1)$。同樣，我們用上划線代表上升冪$n^{\overline{k}}=n(n+1)(n+2) \cdots (n+k-1)$。
斯特林數就是我們轉化的工具，我們有下面兩條公式：
$$ x^n = \sum_{k=0}^{n} \left\{ \begin{matrix}n\\k\end{matrix} \right\} x^{\underline{k}} $$
$$ x^{\overline{n}} = \sum_{k=0}^{n} \left[ \begin{matrix}n\\k\end{matrix} \right] x^k $$
兩條式子都可以利用上面的遞推式通過歸納法證明（證明過程中需要處理下降冪，如$x^{\underline{k+1}}=x^{\underline{k}}(x-k)$，注意和通常冪處理時的差異）。
這兩條式子還有簡單的組合意義，如果不想通過歸納那么暴力地證明的話也可以從這個角度證明。
第一條式子是說，用$x$種顏色為$n$個點染色，方案數顯然是$x^n$；但是我們也可以用另外一個方法表示，那就是枚舉同種顏色的集合是什么（斯特林子集數），然后從$x$個顏色里選出不同的顏色賦給這些集合（下降冪）。
第二條式子是說，用$n$顆珠子串成若干串項鏈，其中每串項鏈上珠子的顏色必須相同，那么我們就可以枚舉項鏈的組成（斯特林輪換數），然后統一給每串項鏈分配顏色（$x^k$）；如果用另外一種方法表示的話，我們考慮按照編號從小到大加入每個珠子，每個珠子可以選擇在$x$種顏色里面選一種並且自己成為一串新的項鏈，也可以選擇接在之前某個珠子后面並繼承一樣的顏色，那么我們的操作方案數是$x^{\overline{n}}$，同時我們發現我們每種操作方案唯一地對應了那些項鏈的組合，因為我們相當於把每條項鏈在編號最小的珠子前剖開，這是唯一的。
第二條式子加以變形可以得到下降冪的版本。這是因為觀察到上升冪和下降冪的式子中剛好是符號反轉，如果我們寫成多項式的話就會發現下降冪除了帶有交錯符號外和上升冪是一樣的。
$$ x^{\underline{n}} = \sum_{k=0}^{n} (-1)^{n-k} \left[ \begin{matrix}n\\k\end{matrix} \right] x^k $$
第一條式子也能變形，不過不是我們的主題。這次我們要用的就是下降冪。
首先我們來解釋一下標題的括號里的內容…我們要做的是離散微積分（有限微積分）。如果了解微積分的話（或者自學），就會知道通過定積分計算函數線下面積這一點其實是一種小技巧，就是小學生都會的裂項求和的原理。通過把$f(x)$寫成$g(x+\Delta x)-g(x)$的形式再把它們加起來，相鄰兩項就會抵消，最后只剩下頭和尾。離散微積分做的也是這種事，只是傳統的微積分里我們取的差是無窮小，而現在我們取$1$。於是我們得到的就是整點處的函數值的和。
為了簡潔易懂，這里就不再引入太多概念了（各種函數的求和、分部求和公式什么的就不說了），我們就說最基本的…
差分算子$\Delta$：我們定義$\Delta f(x) = f(x+1)-f(x)$。
於是對於一個和式$\sum_{k=a}^{b} f(k)$，我們找到一個原函數$g(x)$使得$\Delta g(x)=f(x)$，那么原來的和式就可以寫成：
$$ \sum_{k=a}^{b} g(k+1)-g(k) = g(b+1)-g(a) $$
（如果你原來不了解微積分的話你應該對這種相消表示贊嘆）
所以如果我們能夠方便地求出原函數，就能快速計算所有和式了。不過通常的問題是原函數並不好求，這就是我們選擇下降冪的原因。
因為$ \Delta x^{\underline{n}} = (x+1)^{\underline{n}}-x^{\underline{n}} = [x+1-(x-n+1)]x^{\underline{n-1}} = nx^{\underline{n-1}} $。
這就像通常冪的求導一樣，非常簡單。反過來我們也能輕易地求出原函數就是$\frac{x^{\underline{n+1}}}{n+1}$。
於是我們對我們的自然數冪和問題有想法了，那就是通過斯特林數轉化為下降冪，然后用離散微積分幫助我們求和，然后再轉化回去。我們可以直接寫出式子：
$$ \begin{aligned} \sum_{k=0}^{n} k^t & = \sum_{k=0}^{n} \sum_{i=0}^{t} (-1)^{t-i} \left\{ \begin{matrix}t\\i\end{matrix} \right\} k^{\underline{i}} \\ & = \sum_{i=0}^{t} \left\{ \begin{matrix}t\\i\end{matrix} \right\} \sum_{k=0}^{n} k^{\underline{i}}\\ & = \sum_{i=0}^{t} \left\{ \begin{matrix}t\\i\end{matrix} \right\} \frac{(n+1)^{\underline{i+1}} }{i+1} \\ & = \sum_{i=0}^{t} \left\{ \begin{matrix}t\\i\end{matrix} \right\} \frac{1}{i+1} \sum_{j=0}^{i+1} \left[ \begin{matrix}i+1\\j\end{matrix} \right] (-1)^{i+1-j} (n+1)^j \end{aligned} $$
於是就這樣了，到這個形式不需要繼續化簡就已經可以$O(t^2)$計算了。而且從這里我們可以看出，自然數冪和確實就是一個關於$n$的$t+1$次多項式。同時，如果我們想求得多項式的系數也很容易：注意到這個$n+1$有些扎眼，我們可以把$0$到$n-1$的冪和用這個公式算（那樣這里就是$n$了），然后單獨加上一個$n^t$。而加上$n^t$只是在多項式的$t$次項系數加一而已。所以我們可以$O(t^2)$地得到多項式的系數，比上面那個好多了。
方法3、高階差分(牛頓級數)
在上一部分我們已經定義了一個差分算子啦…高階差分其實就是對一個東西進行多次差分得到一些結果。沒有自己試一試的話當然是覺得很玄乎的，所以不妨拿筆寫一寫三階差分$\Delta ^3 f(x)$是什么。
寫出來就是$\Delta^3 f(x) = f(x+3)-3f(x+2)+3(f+1)-f(x)$，對吧？看到這些系數基本人人都能想到二項式系數吧，如果算了二階和四階差分的話就會更加肯定這個想法。
我們可以得到結論$ \Delta^n f(x) = \sum_{k=0}^{n} \binom{n}{k} (-1)^{n-k} f(x+k) $
如果要嚴謹地證明的話可以用歸納法，或者如果有比較抽代或者線代的思想觀念，可以直接對着算子套用二項式定理得到（比如我們把多項式寫成向量，就能把差分算子寫成矩陣，同時還能寫出另一個讓$f(x)$遞推到$f(x+1)$的矩陣，我們叫做移位算子，我們發現差分矩陣就是移位矩陣減去單位矩陣，於是套用二項式定理就能得到結果）。
這條式子還不能告訴我們什么，接下來的牛頓級數會揭露它的用處：
牛頓級數：$f(x)=c_d \binom{x}{d} + c_{d-1} \binom{x}{d-1} + \cdots + c_1 \binom{x}{1} + c_0 \binom{x}{0}$
我們先來做一些解釋，我們可以唯一地把一個多項式寫成牛頓級數，因為二項式系數可以寫成下降冪除以卷積的形式，而上一部分中斯特林數告訴我們通常冪可以唯一地轉化為下降冪的和，我們把式子重新整理為下降冪的和，然后在每一項系數上乘以一個階乘，於是我們就會得到一個牛頓級數。
牛頓級數的差分具有非常好的性質，因為按照二項式系數的遞推公式可以得到二項式系數的差分$\Delta\left( \binom{x}{k} \right)=\binom{x+1}{k}-\binom{x}{k}=\binom{x}{k-1}$。於是我們經過$i$次差分后，$\binom{x}{0}$的系數就是原來的$c_i$。同時我們還能注意到，當$x$取$0$時，除了這一項以外的其他項全部會變成$0$，所以$\Delta^i f(0)$就等於$c_i$。
現在回過去看上面的式子，它告訴了我們$\Delta^n f(0)$只和$f(0),f(1),\cdots,f(n)$有關，所以如果我們知道了前面這些項的值，我們就可以得到整個多項式了。我們用$O(t^2)$的時間暴力算出差分，然后代回牛頓級數里就能得到原多項式，復雜度還是$O(t^2)$的。至於前$t+2$項的$t$次方和，隨便怎么算都行，總復雜度還是$O(t^2)$。
不過，如果我們的目標只是求值而不是求多項式的系數的話，我們還可以直接把高階差分公式帶進牛頓級數：
$$\begin{aligned} f(n)&=\sum_{k=0}^{d} \binom{n}{k} \sum_{i=0}^{k} \binom{k}{i} (-1)^{k-i} f(i) \\ &= \sum_{k=0}^{d} \sum_{i=0}^{k} \binom{n}{i} \binom{n-i}{k-i} (-1)^{k-i}f(i)\\ &= \sum_{i=0}^{d} f(i)\binom{n}{i} \sum_{k=i}^{d} \binom{k-n-1}{k-i} \\ &= \sum_{i=0}^{d} f(i)\binom{n}{i} \sum_{k=0}^{d-i} \binom{k+i-n-1}{k} \\ &= \sum_{i=0}^{d} f(i)\binom{n}{i} \binom{d-n}{d-i} \\ &= \sum_{i=0}^{d} (-1)^{d-i}f(i)\binom{n}{i} \binom{n-i-1}{d-i} \\ &= \sum_{i=0}^{d} (-1)^{d-i}f(i) \frac{n!}{(n-i)!i!} \cdot \frac{(n-i-1)!}{(d-i)!(n-d-1)!} \\ &= \sum_{i=0}^{d} (-1)^{d-i}f(i) \frac{n^{\underline{d+1}}}{(n-i)i!(d-i)!} \end{aligned}$$
最后得到的那個東西可以$O(d)$算出來，預處理$O(d)$的階乘就可以快速算下面，然后注意到分子其實可以和分母中的$(n-i)$約掉一項，所以分子可以變成兩段連續下降冪的乘積，預處理即可避免除法。對於前$t+2$項的$t$次方和，我們可以使用線性篩，因為冪函數是積性的，對於$O\left(\frac{t}{\ln t}\right)$個質數我們用快速冪算，其它的直接篩出來，這樣總復雜度還是$O(t)$的。於是對於單點求值，我們有了一個非常快的算法。（也有人把這個叫線性插值…不過因為這個詞有別的意思所以我不知道該怎么叫了）
另外注意到我們上面的推導與我們這次要求的自然數冪和沒什么關系，也就是說，這種方法可以非常普遍地適用於多項式的求和。所以還是非常強大的。
方法4、伯努利數(指數生成函數)
這次我們的最后一個方法是專門處理自然數冪和的方法…就是直接用伯努利數。這是由雅各布·伯努利發現的一個數列…
如果我們用別的算法算出對於不同的$t$的自然數冪和的多項式，然后寫在一張紙上，其實可以觀察到一些性質（這里就不列出來了）。雅各布·伯努利就發現了其中的規律…
這里為了方便，我們修改一下我們的記號：$S_t(n) = \sum_{k=0}^{n-1} k^t$。我們把$n^t$的一項去掉了。這會比較方便接下來的分析，並且影響不大，就像我們在斯特林數那一節所做的一樣。
伯努利的結果是：
$$S_t(n) = \frac{1}{t+1} \sum_{k=0}^{t} \binom{t+1}{k}B_k n^{t+1-k}$$
其中$B_k$是一些常數。它們就是伯努利數。這個式子的證明我們先放一放，先看看它有什么用。
注意到，當$n=1$時，$S_t(n)$就變成了$0^t$，所以只有$t=0$時它才是$1$，其它時候都是$0$。注意到在這個時候，我們得到了一個一側含有伯努利數的式子：
$$\sum_{k=0}^{t} \binom{t+1}{k} B_k = 0，t>0$$
而如果$t=0$，我們可以得到$B_0=1$。現在我們就可以遞推算$B_k$了，移項即可得到遞推式，復雜度是$O(t^2)$的。
$$ B_k = -\frac{1}{k+1} \sum_{i=0}^{k-1} \binom{k+1}{i} B_i $$
而且注意到上面伯努利得到的結論直接就是一個多項式的形式，所以通過伯努利數得到多項式是$O(t)$的，是這些方法里最快的了。
接下來我們需要證明這個東西…這個東西可以用歸納法證明，但是看起來非常復雜暴力…好在還有一種方法是指數生成函數。
一個數列$f_n$的指數生成函數是$\hat F(z) = \sum_{n \geq 0} f_n \frac{z^n}{n!} $。
所以如果我們把兩個指數生成函數相乘，得到的指數生成函數是它們的二項卷積的指數生成函數。比如對於$\hat F(z) \hat G(z) = \hat H(z)$：
$$ h_n = \sum_{k=0}^{n} \binom{n}{k} f_k g_{n-k} $$
我們可以對開始那個遞推式使用指數生成函數：我們用$n$代替$t+1$，然后在兩側加上$B_n$得到
$$\sum_{k=0}^{n} \binom{n}{k} B_k = B_n，n>1$$
左側其實是$\hat B(z)$與一個系數全是$1$的數列的指數生成函數($e^z$)二項卷積，同時我們補上$n=1$的情況，此時在右邊會多一個$1$。所以
$$ \hat B(z) e^z = \hat B(z) + z $$
變形就能得到$ \hat B(z) = \frac{z}{e^z-1} $，這就是伯努利數的指數生成函數。（這也得出了一種通過FFT對多項式求逆的一種$O(t \log t)$的預處理伯努利數的方法）
接下來考慮$t$次方和。我們嘗試利用指數生成函數湊$t$次方和。因為普通生成函數下我們得到的是分母上的等差數列，只能寫成調和數的形式，所以實際意義不大。而指數生成函數中我們有希望把等差數列挪到指數上，從而得到一個等比數列，然后就可以化簡了。
$$ e^{kz} = \sum_{i \geq 0} k^i \frac{z^i}{i!} $$
這就是數列$k^0,k^1,k^2,\cdots$的指數生成函數，所以如果我們要求自然數冪和，我們只要把$0$到$n-1$的$e^{kz}$相加：
$$ \sum_{k=0}^{n-1} e^{kz} = \frac{e^{nz}-1}{e^z-1} $$
這就得到了一個簡單的形式了。聯系伯努利數的生成函數就有：
$$ \hat S_t(z) = \hat B(z) \frac{e^{nz}-1}{z} $$
右側又可以變成類似兩個數列的二項卷積的形式，不過有一點不同，其中$\frac{e^{nz}-1}{z}$其實就是把$n^0,n^1,n^2,\cdots$的生成函數去掉常數項再在普通生成函數意義下挪了一位。這就導致了公式的二項式系數里面那個奇怪的$+1$。最后我們取$ \hat S_t(n) $的第$t$項系數再乘以$t!$，就能得到上面那個公式了。
這大概是最針對性的自然數冪和的計算方法了，對於求多項式有很快的速度，超過其它幾種方法，預處理上也有FFT的快速算法，所以效率是相當不錯的。（只是在普適性上，看起來還是上面的高階差分的辦法更加強大一點）