【学习笔记】自然数幂和

温馨提示: 本文文档大小约(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). (“互不相同”指对于任意(0le i<jle n)有(x_i e 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) imes (n+1))的系数矩阵可逆。

实际上系数矩阵就是著名的范德蒙德(Vandermonde)矩阵: (egin{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^nend{bmatrix})

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

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

[egin{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'=egin{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 imes n)的Vandermonde矩阵，由归纳知其行列式为(prod_{1le i<jle n}(a_j-a_i)), 故(|V|=prod_{0le i<jle n}(a_j-a_i)), 证毕。

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

任意多项式的插值

直接给出公式。

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

证明 代入(x_k), 当(i e 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(nlog ^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_iprod_{0le jle n, j e i}frac{n-j}{i-j}=sum^{m}_{i=0}y_ifrac{prod_{0le jle n}(m-j)}{prod_{0le jle n,j e i}(i-j) imes (m-i)}), 可以用阶乘和逆元之类的方法优化，时间复杂度(O(nlog n)-O(n)) 或 (O(n)-O(n)).

第二类斯特林(Stirling)数法

概述

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

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

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

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

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

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

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

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

定理4 (egin{Bmatrix}n\mend{Bmatrix}=frac{1}{m!}sum^{m}_{i=0}(-1)^iegin{pmatrix}m\iend{pmatrix}(m-i)^n)

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

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

分别构造多项式(F(x)=sum_{ige 0}frac{(-1)^i}{i!}x^i)与(G(x)=sum_{ige 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}egin{Bmatrix}n\jend{Bmatrix}i^{underline j}=sum^{n}_{j=0}egin{Bmatrix}n\jend{Bmatrix}sum^{m}_{i=0}i^{underline j}=sum^{n}_{j=0}frac{1}{j+1}egin{Bmatrix}n\jend{Bmatrix}(m+1)^{underline {j+1}} ]

于是我们得到了一种通过第二类斯特林数求自然数幂和的方法。如果(n)不固定，那么可以(O(n^2))预处理第二类斯特林数；如果(n)固定，那么可以(O(nlog n))求第二类斯特林数一整行。求出第二类斯特林数之后，可以在不超过(O(nlog 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(nlog n)).

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

伯努利数(B_n)由以下递归关系定义: (sum^{n}_{i=0}{n+1choose 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_{nge 0}frac{B_n}{n!}=frac{x}{e^x-1}).

证明

[sum^{n}_{i=0}{n+1choose i}B_i=[n=0]\sum^{n+1}_{i=0}{n+1choose i}B_i=B_{n+1}+[n=0]\sum^{n}_{i=0}{nchoose 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(nlog n))时间内预处理伯努利数前(n)项。

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

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

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

[S_m(x)=sum_{nge 0}S(m,n)frac{x^n}{n!}\ =sum_{nge 0}sum^{m-1}_{i=0}i^nfrac{x^n}{n!}\=sum^{m-1}_{i=0}sum_{nge 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_{nge 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(Nlog N))时间内对固定的(m)和每个(n=1,2,...,N)求出(S(m,n)).

三种方法的比较

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

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

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

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

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

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

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