求解自然数幂和的n种方法

问题的引入

给定$n,k$求$$\sum _{i=1}^n{i}^k$$

1. 循环

四年级应该会循环了。

能做到$O(nk)$的优秀时间复杂度。

2. 快速幂

五年级学了快速幂之后就能做到$O(nlo{g}_{2}k)$

请不要小看这个算法。有时候在特定的情况下（例如$n$很小，或$1\to n$的距离变得很小时），这个复杂度真的很优秀。

3. 差分法

六年级应该知道差分和二项式定理了。那么：$$(a+1{)}^k-a^k=\sum _{i=0}^{k-1}{C}_k^i{a}^i$$

于是：
$$\begin{gathered} (n+1{)}^k-1=\sum _{i=1}^n(i+1{)}^k-i^k \\ =\sum _{i=1}^n\sum _{j=0}^{k-1}{C}_k^j{i}^j \\ =\sum _{i=0}^{k-1}{C}_k^i\sum _{j=1}^n{j}^i \\ =\sum _{i=0}^{k-1}{C}_k^{i}S(i) \end{gathered}$$$$\therefore (n+1{)}^{k+1}-1=\sum _{i=0}^k{C}_{k+1}^{i}S(i)$$

把$i=k$时移项，可以得到$$(k+1)S(k)=(n+1{)}^{k+1}-1-\sum _{i=0}^{k-1}{C}_{k+1}^{i}S(i)$$

所以$$S(k)=\frac{(n+1{)}^{k+1}-1-{\sum }_{i=0}^{k-1}{C}_{k+1}^{i}S(i)}{k+1}$$

同时仔细观察这个式子，我们发现，$k$次方和的求和公式是$k+1$次的。归纳证明即可。

3. 倍增

初一应该会倍增了，所以我们令$f_{n,k}={\sum }_{i=1}^n{i}^k$。

当$n$是奇数的时候直接由$f_{n-1,k}+n^k$转移过来。偶数的时候拆开来，运用简单的二项式定理，一波式子推得：$$f(n,k)=f(\frac{n}{2},k)+\sum _{j=0}^k{C}_k^j\ast f(\frac{n}{2},j)\ast {\frac{n}{2}}^{k-j}$$

每一层的$f_{n,k}$我们计算的时间复杂度都是$O(k^2)$的，$lo{g}_n$层，时间复杂度$O(k^{2}lo{g}_n)$.

4. 高斯消元

初一应该会高斯消元了。这是个大脑洞。虽然时间复杂度比上一个还劣一些。

根据$k$次方和的求和公式是$k+1$次的，所以列出$k+2$条式子就可以唯一确定这个多项式。

时间复杂度$O(k^3)$.

5. 第一类斯特林数

初二来学习一下斯特林数。

第一类斯特林数我们一般清楚的是它的组合意义，即把$n$个元素分成$k$个圆排列的方案。根据组合意义，我们不难推出它的式子是$$S_u(n,m)=S_u(n-1,m-1)+(n-1)S_u(n-1,m)$$

但事实上，我们求解自然数幂和需要用到的是它的原始定义：
$$x^{n\downarrow }=x\cdot (x-1)\cdot (x-2)\cdots (x-n+1)=\sum _{k=0}^n{s}_s(n,k)\cdot x^k$$$$x^{n\uparrow }=x\cdot (x+1)\cdot (x+2)\cdots (x+n-1)=\sum _{k=0}^n{s}_u(n,k)\cdot x^k$$

这里需要注意，第一类斯特林数根据定义分成了有符号$S_s$和无符号$S_u$两种。事实上，我们可以很轻松的从这个原始定义推出它的组合意义。

因为$$\sum _{k=0}^n{S}_u(n,k)\cdot x^k=x^{n\uparrow }=x^{(n-1)\uparrow }\cdot (x+n-1)$$$$=\sum _{k=0}^{n-1}{S}_u(n-1,k)x^{k+1}+(n-1)\sum _{k=0}^{n-1}{S}_u(n-1,k)x^k$$

对比两边$x^m$的系数，可以得到$$S_u(n,m)=S_u(n-1,m-1)+(n-1)S_u(n-1,m)$$继续推有符号的，可以得到$$S_s(n,m)=S_s(n-1,m-1)-(n-1)S_s(n-1,m)$$事实上，我们可以完全不用记第一类斯特林数的组合意义，通过公式直接推出来即可。当然记了更好，还可以验证。

所以，根据第一类斯特林数的的定义，得到：$$\prod _{x=0}^{k-1}(n-x)=\sum _{k=0}^n{S}_s(n,k)x^k$$

于是我们可以得到一个显然的式子是：$$n^m=n^{m\downarrow }-\sum _{k=0}^{m-1}{S}_s(m,k)\cdot n^k$$

我们继续推，发现下降幂的和是可以写成一个组合数的形式的，比方说$$\sum _{i=m}^n{i}^{m\downarrow }=\sum _{i=m}^n\frac{i!m!}{(i-m)!m!}=m!\sum _{i=m}^n\left(\genfrac{}{}{0ex}{}{i}{m}\right)=m!\left(\genfrac{}{}{0ex}{}{n+1}{m+1}\right)$$

而后面那一坨式子也是可以化简的，比方说$$\sum _{i=0}^n\sum _{k=0}^{m-1}{S}_s(m,k)\cdot i^k=\sum _{k=0}^{m-1}{S}_s(m,k)\sum _{i=0}^n{i}^k$$

发现后面那条式子${\sum }_{i=0}^n{i}^k$的$k$是降了阶的，所以可以边处理边记录一下，就不用重新算了，时间复杂度就变成了$O(k^2)$。而处理$S_s(m,k)$也是$O(k^2)$级别。

事实上如果当你升入初三，$S_s(m,k)$就可以运用分治$NTT$做到$O(klo{g}^{2}k)$了，虽然然并卵。

但请注意，这种方法虽然时间复杂度是$O(k^2)$级别的，但是它并非没有什么用，因为它——不用做除法

6. 第二类斯特林数

第二类斯特林数的组合意义就是$n$个元素分成$m$个集合，且集合非空的方案数。

基本性质是$$\left\{\begin{array}{c}n\\ m\end{array}\right\}=\left\{\begin{array}{c}n-1\\ m-1\end{array}\right\}+m\cdot \left\{\begin{array}{c}n-1\\ m\end{array}\right\}$$

考虑它的通项公式，可以先把所有集合标号，最后除以集合的阶乘即可，那么考虑容斥，枚举非空集合个数$i$，可以得到$$\left\{\begin{array}{c}n\\ m\end{array}\right\}=\frac{1}{m!}\sum _{i=0}^m(-1{)}^i(\genfrac{}{}{0ex}{}{m}{i})(m-i{)}^n$$

接下来继续推导自然数幂和。

显然！！$$i^k=\sum _{j=0}^k\left\{\begin{array}{c}k\\ j\end{array}\right\}{i}^{j\downarrow }$$

继续推导$$\sum _{i=1}^n{i}^k=\sum _{i=1}^n\sum _{j=0}^k\left\{\begin{array}{c}k\\ j\end{array}\right\}{i}^{j\downarrow }$$$$=\sum _{i=1}^n\sum _{j=0}^k\left\{\begin{array}{c}k\\ j\end{array}\right\}j!\left(\begin{array}{c}i\\ j\end{array}\right)$$$$=\sum _{j=0}^k\left\{\begin{array}{c}k\\ j\end{array}\right\}j!\sum _{i=j}^n\left(\begin{array}{c}i\\ j\end{array}\right)$$$$=\sum _{j=0}^k\left\{\begin{array}{c}k\\ j\end{array}\right\}j!\left(\genfrac{}{}{0ex}{}{n+1}{j+1}\right)$$

很明显，除去预处理第二类斯特林数的复杂度，后面是一样不用做除法的，可以做到$O(k)$.

那么时间复杂度决定于预处理第二类斯特林数的复杂度。显然可以用$O(k^2)$递推。

然而事实上，我们来看看斯特林数的通项公式：$$\left\{\begin{array}{c}n\\ m\end{array}\right\}=\frac{1}{m!}\sum _{i=0}^m(-1{)}^i(\genfrac{}{}{0ex}{}{m}{i})(m-i{)}^n$$

一拼凑，咦~$$\left\{\begin{array}{c}n\\ m\end{array}\right\}=\sum _{i=0}^m\frac{(-1{)}^i}{i!}\cdot \frac{(m-i{)}^n}{(m-i)!}$$

这原来可以写成形如${\sum }_{i=0}^{m}f(i)\ast g(m-i)$的卷积形式。于是第一个多项式的第$i$项系数是$\frac{(-1{)}^i}{i!}$，另一个多项式的第$i$项系数是$\frac{i^n}{i!}$，卷积后第$i$项的系数就是$\left\{\begin{array}{c}n\\ i\end{array}\right\}$.

于是愉快的将时间变成了$O(KlogK)$

7. 差分表

初三来学习一下差分表吧。

对于任何一个序列$a_0,a_1,...,a_n,...$我们都可以定义它的差分序列$\Delta a_0,\Delta a_1,...,\Delta a_n,...$，其中$\Delta a_i=a_{i+1}-a_i$

类似的，我们可以构造序列{Δan}\{\Delta a_n\}{Δan​}的二阶、三阶…$k$阶差分序列。 不妨记为{Δ2an},...,{Δkan}\{\Delta^2 a_n\},...,\{\Delta^k a_n\}{Δ2an​},...,{Δkan​}

令序列是一个$p$次多项式，那么差分表一个很重要且很显然的性质是：$$\forall n>=0,{\Delta }^{p+1}{a}_n=0$$这是由于每次差分都必然会把最高次项消去！

另外一个很重要的性质就是差分表的线性性，即如果$f_n=k_1{g}_n+k_2{h}_n$，那么一定有$$\forall p,n,{\Delta }^p{f}_n=k_1{\Delta }^p{g}_n+k_2{\Delta }^p{h}_n$$

而其最重要的一个性质就是，任何一个$p$阶多项式，都必定可以由其差分表的第一条对角线确定。为了证明这个结论，不妨先考虑最简单的情况：

差分表的一条对角线为$0,...,0,1,0,...$，即第一条对角线上只有第$p$个位置为$1$，其他都为$0$，那么可以写出这个序列的通项公式$$f_n=c(n)(n-1)(n-2)...(n-p+1)$$

代入$n=p,f_p=1$，得到$$c=\frac{1}{p!}$$所以可以得到$$f_n=\frac{n!}{p!(n-p)!}=\left(\genfrac{}{}{0ex}{}{n}{p}\right)$$

那么根据差分表的线性性，我们就可以得知$$f_n=\sum _{i=0}^p{c}_i\left(\genfrac{}{}{0ex}{}{n}{i}\right)$$
由于$$\sum _{k=0}^{n}f(k)=\sum _{k=0}^p{c}_k\left(\genfrac{}{}{0ex}{}{n+1}{k}\right)$$所以利用差分表，我们可以在$O(p^2)$的时间复杂度求解类似于$$\sum _{i=0}^n{f}_i$$的式子。

回到自然数幂和的问题上，我们把$f_i=i^k$代入计算前$p$项的值，通过$p^2$的时间复杂度处理出差分表的第一条对角线，设这个对角线为$c(p,0),c(p,1),c(p,2),\dots ,c(p,p)$，那么答案就是$$\sum _{k=0}^{p}c(p,k)\left(\genfrac{}{}{0ex}{}{n+1}{k}\right)$$

8. 伯努利数

初三再来学一学伯努利数吧。

根据伯努利数的生成函数定义，可知$$\frac{x}{e^x-1}=\sum _{i\ge 0}{B}_i\ast \frac{x^i}{i!}$$

由于$$\frac{x}{e^x-1}\cdot (e^x-1)=x$$

且$$[x^n]\frac{x}{e^x-1}\cdot (e^x-1)=\sum _{i=0}^{n-1}\frac{B_i}{i!}\cdot \frac{1}{(n-i)!}=[n=1]$$

两边都乘上$n!$可以得到$$\sum _{i=0}^{n-1}\left(\genfrac{}{}{0ex}{}{n}{i}\right)B_i=[n=1]$$

这是伯努利数的一个基本性质。后面会用到。

我们再定义一个多项式$B_n(t)$表示$$B_n(t)=\sum _{k=0}^{n-1}{B}_k\ast t^{n-k}\ast C_n^k$$

然后我们发现
$$B_n(t+1)-B_n(t)=\sum _{k=0}^{n-1}{B}_k\ast ⟮(t+1{)}^{n-k}-t^{n-k}⟯C_n^k$$$$=\sum _{k=0}^{n-1}{B}_k\ast (\sum _{i=0}^{n-k-1}{C}_{n-k}^i\ast t^i)\ast C_n^k$$$$=\sum _{k=0}^{n-1}{B}_k\ast (\sum _{i=0}^{n-k-1}{C}_{n-k}^i\ast t^i\ast C_n^k)$$$$=\sum _{i=0}^{n-1}{B}_k\ast \sum _{i=0}^{n-k-1}\frac{n!t^i}{i!k!(n-k-i)!}$$$$=\sum _{k=0}^{n-1}{B}_k\ast (\sum _{i=0}^{n-k-1}{C}_{n-i}^k\ast t^i\ast C_n^i)$$$$=\sum _{i=0}^{n-1}{C}_n^i\ast t^i\ast \sum _{k=0}^{n-1-i}{B}_k\ast C_{n-i}^k$$

注意到后面只有当$n-i-1=0$时值为$1$，于是$$B_n(t+1)-B_n(t)=n\ast t^{n-1}$$

然后我们考虑差分，就有$$\sum _{t=0}^{n-1}{B}_k(t+1)-B_k(t)=k\cdot \sum _{i=0}^{n-1}{i}^{k-1}$$即$$B_{k+1}(n+1)=(k+1)\cdot \sum _{i=0}^n{i}^k$$

可得自然数幂和$$\sum _{i=0}^n{i}^k=\frac{1}{k+1}\cdot \sum _{i=0}^k{B}_i{n}^{k+1-i}\left(\genfrac{}{}{0ex}{}{k+1}{i}\right)$$

问题转化成了求$B_i$，注意到它的生成函数定义，事实上我们只需要求$$\sum _{i\ge 0}\frac{x^i}{(i+1)!}$$在模$x^{k+1}$的逆元即可。

时间复杂度$O(KLogK)$，当然如果递推的话，也是可以轻松做到$O(k^2)$的。

9. 拉格朗日插值法

两大作用：

1. 快速根据点值逼近原函数.

2. 取点值对大于$n$时唯一确定$n$次多项式.

Example

例如对于$$\sum _{i=1}^{n}i=\frac{n(n+1)}{2}$$

我们知道它的通项公式是二次的，所以我们只需要三个点值对就可以唯一确定这个多项式：$$(1,1),(2,3),(3,6)$$

General method

对于已知的$n+1$个点对$(x_0,y_0),(x_1,y_1)...(x_n,y_n)$，求$n+1$个函数$f_i$，使得该函数在$x_i$处取得对应的$y_i$值，其余$x_j$处为$0$，最后把这$n+1$个函数线性结合即可。

$$f_i(x)=\frac{\prod _{j̸=i}(x-x_j)}{\prod _{j̸=i}(x_i-x_j)}\ast y_i$$$$g(x)=\sum _{i=0}^n{f}_i(x)$$

Practice

例如我们要求自然数幂和.

各种方法可以证明$i^k$的和是$k+1$次的, 所以我们只需要给出$k+2$个点值表达,就可以求得通项公式.

记$S(n)={\sum }_{i=1}^n{i}^k$, 则
$$S(n)=\sum _{i=1}^{k+2}{y}_i\prod _{j=1,i̸=j}^{k+2}\frac{n-x_j}{x_i-x_j}=\sum _{i=1}^{k+2}{y}_i\frac{{\prod }_{j=1,i̸=j}^{k+2}(n-j)}{{\prod }_{j=1,i̸=j}^{k+2}(i-j)}$$

那么时间复杂度就在预处理$y_i$上面了, 利用线筛,可以做到$O(k)$级别.

后面的那一部分可以预处理,具体的说,就有:$$S(n)=\sum _{i=1}^{k+2}{y}_{i}pre[i-1]suf[i+1][(-1{)}^{k+2-i}(i-1)!(k+2-i)!{]}^{-1}$$

$Code$

时间复杂度：$O(\frac{k}{l{n}_k}lo{g}_{2}k)=O(k)$.

另外，注意逆元要预处理，实测：$k=1e7$时$0.9s$，可以说是非常优秀了

#include <bits/stdc++.h> 

#define F(i,a,b) for (int i = a; i <= b; i ++)
#define G(i,a,b) for (int i = a; i >= b; i --)

const int Mo = 998244353, M = 1e6 + 10;

using namespace std;

int l, r, k, m, y[M], z[M], jc[M], suf[M], pre[M];
bool bz[M];

int ksm(int x, int y) {
	int ans = 1;
	for (; y; y >>= 1, x = (1ll * x * x) % Mo)
		if (y & 1)
			ans = (1ll * ans * x) % Mo;
	return ans;
}

void Init() {
	scanf("%d%d%d", &l,&r,&k), y[1] = 1, m = k + 2;
	F(i, 2, m) {
		if (!bz[i])
			z[++ z[0]] = i, y[i] = ksm(i, k);
		F(j, 1, z[0]) {
			if (z[j] * i > m) break;
			bz[z[j] * i] = 1;
			y[z[j] * i] = (1ll * y[z[j]] * y[i]) % Mo;
			if (i % z[j] == 0) break;
		}
	}
	F(i, 2, m)
		y[i] = (y[i - 1] + y[i]) % Mo;
	jc[0] = 1;
	F(i, 1, m)
		jc[i] = 1ll * jc[i - 1] * i % Mo;
	jc[m] = ksm(jc[m], Mo - 2);
	G(i, m - 1, 1)
		jc[i] = 1ll * jc[i + 1] * (i + 1) % Mo;
}

int Solve(int n) {
	pre[0] = suf[m + 1] = 1;
	F(i, 1, m)
		pre[i] = 1ll * pre[i - 1] * (n - i) % Mo;
	G(i, m, 1)
		suf[i] = 1ll * suf[i + 1] * (n - i) % Mo;

	int Ans = 0;
	F(i, 1, m)
		Ans = (Ans + 1ll * y[i] * pre[i - 1] % Mo * suf[i + 1] % Mo * (((k-i+2)&1) ? (-1) : 1) * jc[i - 1] % Mo * jc[k + 2 - i] % Mo) % Mo;
	return Ans;
}

int main() {
	Init();

	printf("%d\n", (Solve(r) - Solve(l - 1) + Mo) % Mo);
}