温馨提示: 本文文档大小约\(11KB\).

引入

自然数幂和是一个我们从小就耳熟能详的经典问题。定义\(S(m,n)=\sum^{m}_{i=0} i^n\), 显然\(S(m,n)\)为关于\(m\)的不超过\((n+1)\)次多项式，那么给定\(n\), 如何快速求出这个多项式的系数？或者给定\(m\)和\(n\), 如何快速求出\(S(m,n)\)? 这就是本文的讨论内容。

本文介绍三种最常用的方法——拉格朗日插值法、第二类斯特林数法和伯努利数法，这三种方法在解决不同的问题上各有优劣。
当然，还有一些其他的方法，例如组合数递推法、差分斯特林数法等等，本文略过。

符号约定

\(x^{\underline n}\)表示\(x\)的\(n\)次下降幂，即\(x^{\underline n}=\prod^{n}_{i=1}(x-i+1)\).

本文中提到的所有\(i\), 与虚数均无关系。

本文中所有运用数学归纳法进行证明的定理，归纳奠基均省略，请读者自行处理。

拉格朗日(Lagrange)插值法

概述

Lagrange插值法的应用范围很广，给定任意\((n+1)\)个互不相同的点值\((x_i,y_i)\) (\(i=0,1,...,n\))，其可以在\(O(n^2)\)时间内求出一个不超过\(n\)次的多项式\(f(x)\)，满足对于任意\(i\), \(f(x_i)=y_i\). (“互不相同”指对于任意\(0\le i<j\le n\)有\(x_i\ne x_j\))

暴力做法

设\(f(x_i)=\sum^{n}_{i=0} a_ix^i\). 根据\((n+1)\)个点值，我们可以列出\((n+1)\)个方程，第\(j\)个形如\(\sum^{n}_{i=0}x_j^ia_i=y_j\).

直接高斯消元求解，时间复杂度\(O(n^3)\).

插值多项式的唯一性

定理1 给定任意\((n+1)\)对互不相同的点值\((x_i,y_i)\)，恰好存在一个不超过\(n\)次的多项式\(f(x)\)，满足对于任意\(i\), \(f(x_i)=y_i\).

证明 根据暴力做法，我们就是要证明方程左边的\((n+1)\times (n+1)\)的系数矩阵可逆。

实际上系数矩阵就是著名的范德蒙德(Vandermonde)矩阵: \(\begin{bmatrix} 1&1&1&...&1\\ a_0&a_1&a_2&...&a_n\\ a_0^2&a_1^2&a_2^2&...&a_n^2\\ & & & ... & \\ a_0^n&a_1^n&a_2^n&...&a_n^n\end{bmatrix}\)

可以证明，该矩阵行列式为\(\prod_{0\le i<j\le n} (a_j-a_i)\).

考虑数学归纳法，自下而上每一行减去上一行的\(a_0\)倍，该矩阵变为

\[\begin{bmatrix}1&1&1&...&1\\ 0&a_1-a_0&a_2-a_0&...&a_n-a_0\\ 0&a_1(a_1-a_0)&a_2(a_2-a_0)&...&a_n(a_n-a_0)\\&&&...&\\0&a_1^{n-1}(a_1-a_0)&a_2^{n-1}(a_2-a_0)&...&a_n^{n-1}(a_n-a_0)\end{bmatrix} \]

然后对第一列展开，每一行提取公因式之后可得\(|V|=|V'|\prod^{n}_{i=1}(a_i-a_0)\), 其中

\[V'=\begin{bmatrix}1&1&...&1\\a_1&a_2&...&a_n\\a_1^2&a_2^2&...&a_n^2\\&&...&\\a_1^{n-1}&a_2^{n-1}&...&a_n^{n-1}\end{bmatrix} \]

是一个\(n\times n\)的Vandermonde矩阵，由归纳知其行列式为\(\prod_{1\le i<j\le n}(a_j-a_i)\), 故\(|V|=\prod_{0\le i<j\le n}(a_j-a_i)\), 证毕。

由上可知，由于任意\(1\le i<j\le n\)满足\(a_i\ne a_j\), 故该矩阵行列式不为\(0\).

任意多项式的插值

直接给出公式。

定理2 符合\((n+1)\)个点值\((x_i,y_i)\)的多项式为\(f(x)=\sum^{n}_{i=0} y_i\prod_{0\le j\le n,j\ne i}\frac{x-x_j}{x_i-x_j}\).

证明 代入\(x_k\), 当\(i\ne k\)时后面的乘法有一项因数是\(x-x_k=0\), 故该积恒为\(0\). 当\(i=k\)时显然值为\(y_i\). 因此该多项式符合\((n+1)\)个给定点值。

至于这个公式具体是怎么推出来的，大概是考虑求Vandermonde矩阵的逆矩阵，使用矩阵分块法解决，此处不再赘述。

有了这个公式，我们可以先预处理出\(\prod^n_{i=0}(x-x_i)\)的每一项系数(这是一个\((n+1)\)次多项式)，然后对于每个\(i\)暴力用这个多项式除以\((x-x_i)\). 由于多项式除以一次式可在\(O(n)\)时间内完成，该算法时间复杂度为\(O(n^2)\).

使用FFT最好可以做到\(O(n\log ^2n)\), 但是我不会。

那么对于自然数幂和问题，我们可以求出\(S(m,n)\)在\(0,1,2,...,n\)处的点值，然后利用Lagrange插值求出这个\((n+1)\)次多项式的每一项系数，时间复杂度\(O(n^2)\). 如果给定一个\(m\)要求\(S(m,n)\)的值，那么可以直接把\(n\)代入多项式求值，时间复杂度\(O(n^2)-O(n)\) (分别表示预处理和单次询问复杂度), 与\(m\)无关。

对于自然数幂和问题的优化

如果给定\(m,n\)要求\(S(m,n)\), 不需要求出多项式每一项的系数，是否可以做到更优的复杂度？(假设\(n>k\))

\(f(n)=\sum^{m}_{i=0}y_i\prod_{0\le j\le n, j\ne i}\frac{n-j}{i-j}=\sum^{m}_{i=0}y_i\frac{\prod_{0\le j\le n}(m-j)}{\prod_{0\le j\le n,j\ne i}(i-j)\times (m-i)}\), 可以用阶乘和逆元之类的方法优化，时间复杂度\(O(n\log n)-O(n)\) 或 \(O(n)-O(n)\).

第二类斯特林(Stirling)数法

概述

第二类斯特林数是一类重要的计数序列，又称“斯特林子集数”，其定义为\(\begin{Bmatrix}n\\k\end{Bmatrix}\)表示将\(n\)个有编号的物品放入\(k\)个无编号的集合中，每个集合都非空的方案数。

由定义可以得出其递推式: \(\begin{Bmatrix}n\\k\end{Bmatrix}=k\begin{Bmatrix}n-1\\k\end{Bmatrix}+\begin{Bmatrix}n-1\\k-1\end{Bmatrix}\). 考虑最后一个物品，如果将其单独放入一个集合中，则方案数为\(\begin{Bmatrix}n-1\\k-1\end{Bmatrix}\), 否则所有集合都已经有元素了，所以要加以区分，放入不同的集合算不同的方案，方案数为\(k\begin{Bmatrix} n-1\\k\end{Bmatrix}\).

下面两个重要公式给出了第二类Stirling数与幂运算之间的关系。

定理3 \(m^n=\sum^{n}_{i=0} \begin{Bmatrix}n\\i\end{Bmatrix}m^{\underline i}\).

证明 我们可以从代数和组合等不同角度来证明。

证法一 代数做法: 考虑使用数学归纳法。对\(n\)施归纳，

\[m^{n+1}=\sum^{n}_{i=0}m\begin{Bmatrix}n\\i\end{Bmatrix}m^{\underline i}\\ =\sum^{n}_{i=0}(m-i)\begin{Bmatrix}n\\i\end{Bmatrix}m^{\underline i}+\sum^{n}_{i=0}i\begin{Bmatrix}n\\i\end{Bmatrix}m^{\underline i}\\ =\sum^{n}_{i=0}\begin{Bmatrix}n\\i\end{Bmatrix}m^{i+1}+\sum^{n}_{i=0}i\begin{Bmatrix}n\\i\end{Bmatrix}m^{\underline i}\\ =\sum^{n+1}_{i=1}\begin{Bmatrix}n\\i-1\end{Bmatrix}m^{\underline i}+\sum^{n}_{i=0}i\begin{Bmatrix}n\\i\end{Bmatrix}m^{\underline i}\\ =\sum^{n+1}_{i=0}\begin{Bmatrix}n+1\\i\end{Bmatrix}m^{\underline i} \]

证法二 组合做法: \(m^n\)相当于把\(n\)个有标号数放入\(m\)个有标号集合的方案数。这些集合有的是空集，有的是非空集合，考虑枚举有多少个非空集合，先从有标号集合里有顺序地选出这\(i\)个非空集合(方案数\(m^{\underline i}\)), 再计算\(n\)个数放进去的方案数。

定理4 \(\begin{Bmatrix}n\\m\end{Bmatrix}=\frac{1}{m!}\sum^{m}_{i=0}(-1)^i\begin{pmatrix}m\\i\end{pmatrix}(m-i)^n\)

证明 组合意义上相当于上式套容斥，代数意义上相当于上式套二项式反演。

根据定理4，我们发现固定\(n\)之后可以用FFT在\(O(n\log n)\)时间内对每个\(m\)求出第二类斯特林数。

分别构造多项式\(F(x)=\sum_{i\ge 0}\frac{(-1)^i}{i!}x^i\)与\(G(x)=\sum_{i\ge 0}\frac{i^n}{i!}x^i\), 卷积相乘即可。

利用第二类斯特林数求自然数幂和

定理4指出一种快速求第二类斯特林数的方法，而定理3指出第二类斯特林数与自然数的幂的联系。

既然定理3的式子是把通常幂转化为下降幂，那么求自然数下降幂的和是否容易一些呢？

这是一个小学数学问题。

\[\sum^{m}_{i=0}i^{\underline n}=\frac{1}{n+1}\sum^{m}_{i=0}(i+1)^{\underline {n+1}}-i^\underline {n+1}=\frac{(m+1)^{\underline {n+1}}-0^{\underline {n+1}}}{n+1}=\frac{1}{n+1}(m+1)^{\underline {n+1}} \]

那么考虑将自然数幂和转化为自然数下降幂和:

\[S(m,n)=\sum^{m}_{i=0}i^n=\sum^{m}_{i=0}\sum^{n}_{j=0}\begin{Bmatrix}n\\j\end{Bmatrix}i^{\underline j}=\sum^{n}_{j=0}\begin{Bmatrix}n\\j\end{Bmatrix}\sum^{m}_{i=0}i^{\underline j}=\sum^{n}_{j=0}\frac{1}{j+1}\begin{Bmatrix}n\\j\end{Bmatrix}(m+1)^{\underline {j+1}} \]

于是我们得到了一种通过第二类斯特林数求自然数幂和的方法。如果\(n\)不固定，那么可以\(O(n^2)\)预处理第二类斯特林数；如果\(n\)固定，那么可以\(O(n\log n)\)求第二类斯特林数一整行。求出第二类斯特林数之后，可以在不超过\(O(n\log P)\)的时间复杂度内求出\(S(m,n)\), 其中\(P\)为模数，\(O(\log P)\)为求乘法逆元的复杂度。

特别地，如果模数性质不好，没法求逆元怎么办？此时Lagrange插值法失去了用武之地，而该方法依然可用，但时间复杂度退化为\(O(n^2)\). 我们只需要求\((j+1)\)的逆元，而\((m+1)^{\underline {j+1}}\)是把\((j+1)\)个连续整数相乘，肯定有一个是\((j+1)\)的倍数，每次找到那个倍数即可。

伯努利(Bernoulli)数法

注意在这一部分中，我们需要改变\(S(m,n)\)的定义。现在\(S(m,n)\)定义为\(\sum^{m-1}_{i=0}i^n\), 即比原来少了\(m^n\).

概述

伯努利数法与上面两种方法的适用对象有所不同。伯努利数法适用于对于某个固定的\(m\), 以及\(1\)到\(N\)内的所有的\(n\), 求出\(S(m,n)\)的值, 时间复杂度可以做到\(O(n\log n)\).

伯努利数本身就是由伯努利观察自然数幂和的过程中发现的。他对于较小的\(n\)求出了\(S(m,n)\)多项式的每一次项的系数，然后发现\((n+1)\)次系数总是\(1\), 对于\(0\le k\le n\), \((n-k)\)次项系数总是等于一个常数乘以\(n\)的一个下降幂。于是他定义了伯努利数。

伯努利数\(B_n\)由以下递归关系定义: \(\sum^{n}_{i=0}{n+1\choose i}B_i=[n=0]\). 伯努利数的前\(5\)项是: \(B_0=0, B_1=-\frac{1}{2},B_2=\frac{1}{6},B_3=0,B_4=-\frac{1}{30},B_5=0\).

定理5 伯努利数的指数生成函数(EGF)为\(B(x)=\sum_{n\ge 0}\frac{B_n}{n!}=\frac{x}{e^x-1}\).

证明

\[\sum^{n}_{i=0}{n+1\choose i}B_i=[n=0]\\\sum^{n+1}_{i=0}{n+1\choose i}B_i=B_{n+1}+[n=0]\\\sum^{n}_{i=0}{n\choose i}B_i=B_n+[n=1]\\\frac{B_n}{n!}+[n=1]=\sum^{n}_{i=0}\frac{1}{(n-i)!}\frac{B_i}{i!}\\B(x)+x=B(x)e^x\\B(x)=\frac{x}{e^x-1} \]

由此，用FFT多项式求逆就可以在\(O(n\log n)\)时间内预处理伯努利数前\(n\)项。

利用伯努利数求自然数幂和

定理6 伯努利数满足关系式\(S(m,n)=\frac{1}{n+1}\sum^{n}_{i=0}{n+1\choose i}B_im^{n+1-i}\).

证明 考虑固定\(m\), 写出\(S(m,n)\)的指数生成函数\(S_m(x)\).

\[S_m(x)=\sum_{n\ge 0}S(m,n)\frac{x^n}{n!}\\ =\sum_{n\ge 0}\sum^{m-1}_{i=0}i^n\frac{x^n}{n!}\\=\sum^{m-1}_{i=0}\sum_{n\ge 0}\frac{i^n}{n!}x^n\\=\sum^{m-1}_{i=0}e^{ix}\\=\frac{e^{mx}-1}{e^x-1}\\=B(x)\frac{e^{mx}-1}{x}\\=B(x)\sum_{n\ge 0}\frac{m^{n+1}}{(n+1)!}x^n\\ [x^n]S_m(x)=\sum^{n}_{i=0}\frac{B_i}{i!}\frac{m^{n+1-i}}{(n+1-i)!} \]

然后由于是指数生成函数所以\(S(m,n)=n![x^n]S_m(x)\), 证毕。

那么，我们用与上面几例类似的方法，可以用FFT在\(O(N\log N)\)时间内对固定的\(m\)和每个\(n=1,2,...,N\)求出\(S(m,n)\).

三种方法的比较

三种方法适用的问题形式不同，各自的效率也有优劣。

当数据范围可以接受\(O(n^2)\)的时间复杂度时，三种做法都可以不用FFT实现，代码难度都较小。

求出多项式的每一次项系数，这是Lagrange插值法的专长。

当固定\(n\)不固定\(m\)时，Lagrange插值法和第二类Stirling数法较为适用。

当固定\(m\)不固定\(n\)时，Bernoulli数法较为适用。

当模数性质不好无法求逆元时，Stirling数法是较好的选择。

解决有关具体问题时，一定要注意具体的限制(比如询问形式、询问次数、\(n,m\)的大小、是否固定)，选择较好的方法。