本篇口胡写给我自己这样的东西都忘光的残废选手 以及暂时还不会自然数幂和的人…
这里大概给出最简单的几种方法：扰动法(化为递推式)，斯特林数(离散微积分)，高阶差分(牛顿级数)，伯努利数(指数生成函数)…
不同方法的思维难度、普适程度、实现难度、时间复杂度上面都有差异…同时自然数幂和是探究各种求和方法的经典例子，了解多一点它的做法对于处理各种求和问题是有所帮助的…
问题：求$sum_{k=0}^{n} k^t$，其中$t in mathbb{N}$是一个常数。要求求解的时间复杂度与$n$无关。（特别地，$0^0=1$）
方法0、 暴力
我们for一遍$0$到$n$，然后for一遍$0$到$t$就能算出来啦！比较慢？我们可以用快速幂算$k^t$，就不用for一次啦！时间复杂度$Oleft( n log t ight)$。
以上纯属卖萌。
这个复杂度是不能接受的，因为它与$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} inom{t}{i} k^i $。
已经很类似了，接着我们交互求和次序：
$$egin{aligned} S_t(n)+(n+1)^t & = sum_{i=0}^{t} inom{t}{i} sum_{k=0}^{n} k^i \ & = sum_{i=0}^{t} inom{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} = inom{t+1}{t} S_t(n) + sum_{i=0}^{t-1} inom{t+1}{i} S_i(n) $$
移项得到：
$$ S_t(n) = frac{1}{t+1} left( (n+1)^{t+1} - sum_{i=0}^{t-1} inom{t+1}{i} S_i(n) ight) $$
这便是一个我们可以利用的递推式了。我们知道$S_0(n)$就是$0$到$n$的自然数个数$n+1$，我们甚至可以不用知道$S_1(n)$的结论（虽然我们已经知道那是$n(n+1)/2$），因为我们可以递推得到。所以对于一个确定的$n$，我们从$S_0(n)$开始一直递推到$S_t(n)$就可以了。这样子我们需要的时间是$Oleft(t^2 ight)$的，因为每项递推都需要用到前面所有项的信息。这其实挺不错了，因为它与$n$无关，而且也不是太大。
在后面的推导中我们可以发现自然数幂和其实是一个关于$n$的$t+1$次多项式。我们这个做法对于求多项式的系数其实不是很优的，需要$Oleft(t^3 ight)$的时间（因为不再是代一个数算而是代一个多项式算了）。如果知道这是一个多项式的话暴力高斯消元求系数甚至都可以做到$Oleft( t^3 ight)$。所以我们来看下面这个做法。
方法2、斯特林数(离散微积分)
你问那个括号里是什么意思？其实就是对于方法的补充说明啦…
虽然一开始看起来比较高端，但是这个方法是所有方法里面最无脑的了…非常简单，只要知道结论就会了…
首先简单说一下斯特林数…斯特林数有两类，第一类是轮换数，代表$n$个数排成$k$个非空的环的方案数，我们用$left[ egin{matrix}n\kend{matrix} ight]$表示；第二类是子集数，代表$n$个数形成$k$个非空集合的方案数，我们用$left{ egin{matrix}n\kend{matrix} ight}$表示。
根据意义我们可以得到它们的递推式和边界条件。显然如果$n<k$，它们都是$0$；如果$n=k$，它们都是$1$；如果$n>0$且$k=0$，它们都是$0$。
加入一个新数的时候，可以考虑自己形成一个新的环，或者插入任意一个已有的数后面，所以
$$left[ egin{matrix}n\kend{matrix} ight]=left[ egin{matrix}n-1\k-1end{matrix} ight]+(n-1)left[ egin{matrix}n-1\kend{matrix} ight]$$
同样，加入一个新数的时候，可以考虑自己形成一个新的集合，或者插入任意一个已有的集合。
$$left{ egin{matrix}n\kend{matrix} ight}=left{ egin{matrix}n-1\k-1end{matrix} ight}+kleft{ egin{matrix}n-1\kend{matrix} ight}$$
我这次要用它们在我们使用的通常幂$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{ egin{matrix}n\kend{matrix} ight} x^{underline{k}} $$
$$ x^{overline{n}} = sum_{k=0}^{n} left[ egin{matrix}n\kend{matrix} ight] 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[ egin{matrix}n\kend{matrix} ight] 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}$。
于是我们对我们的自然数幂和问题有想法了，那就是通过斯特林数转化为下降幂，然后用离散微积分帮助我们求和，然后再转化回去。我们可以直接写出式子：
$$ egin{aligned} sum_{k=0}^{n} k^t & = sum_{k=0}^{n} sum_{i=0}^{t} (-1)^{t-i} left{ egin{matrix}t\iend{matrix} ight} k^{underline{i}} \ & = sum_{i=0}^{t} left{ egin{matrix}t\iend{matrix} ight} sum_{k=0}^{n} k^{underline{i}}\ & = sum_{i=0}^{t} left{ egin{matrix}t\iend{matrix} ight} frac{(n+1)^{underline{i+1}} }{i+1} \ & = sum_{i=0}^{t} left{ egin{matrix}t\iend{matrix} ight} frac{1}{i+1} sum_{j=0}^{i+1} left[ egin{matrix}i+1\jend{matrix} ight] (-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} inom{n}{k} (-1)^{n-k} f(x+k) $
如果要严谨地证明的话可以用归纳法，或者如果有比较抽代或者线代的思想观念，可以直接对着算子套用二项式定理得到（比如我们把多项式写成向量，就能把差分算子写成矩阵，同时还能写出另一个让$f(x)$递推到$f(x+1)$的矩阵，我们叫做移位算子，我们发现差分矩阵就是移位矩阵减去单位矩阵，于是套用二项式定理就能得到结果）。
这条式子还不能告诉我们什么，接下来的牛顿级数会揭露它的用处：
牛顿级数：$f(x)=c_d inom{x}{d} + c_{d-1} inom{x}{d-1} + cdots + c_1 inom{x}{1} + c_0 inom{x}{0}$
我们先来做一些解释，我们可以唯一地把一个多项式写成牛顿级数，因为二项式系数可以写成下降幂除以卷积的形式，而上一部分中斯特林数告诉我们通常幂可以唯一地转化为下降幂的和，我们把式子重新整理为下降幂的和，然后在每一项系数上乘以一个阶乘，于是我们就会得到一个牛顿级数。
牛顿级数的差分具有非常好的性质，因为按照二项式系数的递推公式可以得到二项式系数的差分$Deltaleft( inom{x}{k} ight)=inom{x+1}{k}-inom{x}{k}=inom{x}{k-1}$。于是我们经过$i$次差分后，$inom{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)$。
不过，如果我们的目标只是求值而不是求多项式的系数的话，我们还可以直接把高阶差分公式带进牛顿级数：
$$egin{aligned} f(n)&=sum_{k=0}^{d} inom{n}{k} sum_{i=0}^{k} inom{k}{i} (-1)^{k-i} f(i) \ &= sum_{k=0}^{d} sum_{i=0}^{k} inom{n}{i} inom{n-i}{k-i} (-1)^{k-i}f(i)\ &= sum_{i=0}^{d} f(i)inom{n}{i} sum_{k=i}^{d} inom{k-n-1}{k-i} \ &= sum_{i=0}^{d} f(i)inom{n}{i} sum_{k=0}^{d-i} inom{k+i-n-1}{k} \ &= sum_{i=0}^{d} f(i)inom{n}{i} inom{d-n}{d-i} \ &= sum_{i=0}^{d} (-1)^{d-i}f(i)inom{n}{i} inom{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$次方和，我们可以使用线性筛，因为幂函数是积性的，对于$Oleft(frac{t}{ln t} ight)$个质数我们用快速幂算，其它的直接筛出来，这样总复杂度还是$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} inom{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} inom{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} inom{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} inom{n}{k} f_k g_{n-k} $$
我们可以对开始那个递推式使用指数生成函数：我们用$n$代替$t+1$，然后在两侧加上$B_n$得到
$$sum_{k=0}^{n} inom{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的快速算法，所以效率是相当不错的。（只是在普适性上，看起来还是上面的高阶差分的办法更加强大一点）