求解自然数幂和的若干种方法
问题的引入
给定\(n,k\)求$$\sum_{i=1}nik$$
1. 循环
四年级应该会循环了。
能做到\(O(nk)\)的优秀时间复杂度。
2. 快速幂
五年级学了快速幂之后就能做到\(O(nlog_2k)\)
请不要小看这个算法。有时候在特定的情况下（例如\(n\)很小，或\(1\rightarrow n\)的距离变得很小时），这个复杂度真的很优秀。
3. 差分法
六年级应该知道差分和二项式定理了。那么：$$(a+1)k-ak=\sum_{i=0}{k-1}C_kia^i$$
于是：
把\(i=k\)时移项，可以得到$$(k+1)S(k)=(n+1){k+1}-1-\sum_{i=0}C_{k+1}^iS(i)$$
所以$$S(k)=\frac{(n+1){k+1}-1-\sum_{i=0}C_{k+1}^iS(i)}{k+1}$$
同时仔细观察这个式子，我们发现，\(k\)次方和的求和公式是\(k+1\)次的。归纳证明即可。
3. 倍增
初一应该会倍增了，所以我们令\(f_{n,k}=\sum_{i=1}^ni^k\)。
当\(n\)是奇数的时候直接由\(f_{n-1,k}+n^k\)转移过来。偶数的时候拆开来，运用简单的二项式定理，一波式子推得：$$f(n,k)=f(\frac{n}{2},k)+\sum_{j=0}kC_kjf(\frac{n}{2},j)\frac{n}{2}^{k-j}$$
每一层的\(f_{n,k}\)我们计算的时间复杂度都是\(O(k^2)\)的，\(log_n\)层，时间复杂度\(O(k^2log_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)$$
但事实上，我们求解自然数幂和需要用到的是它的原始定义：
这里需要注意，第一类斯特林数根据定义分成了有符号\(S_s\)和无符号\(S_u\)两种。事实上，我们可以很轻松的从这个原始定义推出它的组合意义。
因为$$\sum_{k=0}nS_u(n,k)·xk=x{n\uparrow}=x·(x+n-1)$$$$=\sum_{k=0}{n-1}S_u(n-1,k)x+(n-1)\sum_{k=0}{n-1}S_u(n-1,k)xk$$
对比两边\(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}nS_s(n,k)x^k$$
于是我们可以得到一个显然的式子是：$$nm=n-\sum_{k=0}{m-1}S_s(m,k)·nk$$
我们继续推，发现下降幂的和是可以写成一个组合数的形式的，比方说$$\sum_{i=m}ni=\sum_{i=m}n\frac{i!m!}{(i-m)!m!}=m!\sum_{i=m}n\binom{i}{m}=m!\binom{n+1}{m+1}$$
而后面那一坨式子也是可以化简的，比方说$$\sum_{i=0}n\sum_{k=0}S_s(m,k)·ik=\sum_{k=0}S_s(m,k)\sum_{i=0}nik$$
发现后面那条式子\(\sum_{i=0}^ni^k\)的\(k\)是降了阶的，所以可以边处理边记录一下，就不用重新算了，时间复杂度就变成了\(O(k^2)\)。而处理\(S_s(m,k)\)也是\(O(k^2)\)级别。
事实上如果当你升入初三，\(S_s(m,k)\)就可以运用分治\(NTT\)做到\(O(klog^2k)\)了，虽然然并卵。
但请注意，这种方法虽然时间复杂度是\(O(k^2)\)级别的，但是它并非没有什么用，因为它——不用做除法
6. 第二类斯特林数
第二类斯特林数的组合意义就是\(n\)个元素分成\(m\)个集合，且集合非空的方案数。
基本性质是$$\begin{Bmatrix}n\m\end{Bmatrix}=\begin{Bmatrix}n-1\m-1\end{Bmatrix}+m\cdot \begin{Bmatrix}n-1\m\end{Bmatrix}$$
考虑它的通项公式，可以先把所有集合标号，最后除以集合的阶乘即可，那么考虑容斥，枚举非空集合个数\(i\)，可以得到$$\begin{Bmatrix}n\m\end{Bmatrix}=\frac 1 {m!}\sum_{i=0}m(-1)i\binom mi(m-i)^n$$
接下来继续推导自然数幂和。
显然！！$$ik=\sum_{j=0}k\left{\begin{array}{c}{k}\{j}\end{array}\right}i^{j\downarrow}$$
继续推导$$\sum_{i=1}nik=\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!\binom{n+1}{j+1}$$
很明显，除去预处理第二类斯特林数的复杂度，后面是一样不用做除法的，可以做到\(O(k)\).
那么时间复杂度决定于预处理第二类斯特林数的复杂度。显然可以用\(O(k^2)\)递推。
然而事实上，我们来看看斯特林数的通项公式：$$\begin{Bmatrix}n\m\end{Bmatrix}=\frac 1 {m!}\sum_{i=0}m(-1)i\binom mi(m-i)^n$$
一拼凑，咦~$$\begin{Bmatrix}n\m\end{Bmatrix}=\sum_{i=0}m\frac{(-1)i}{i!}\cdot\frac{(m-i)^n}{(m-i)!}$$
这原来可以写成形如\(\sum_{i=0}^mf(i)*g(m-i)\)的卷积形式。于是第一个多项式的第\(i\)项系数是\(\frac {(-1)^i}{i!}\)，另一个多项式的第\(i\)项系数是\(\frac {i^n}{i!}\)，卷积后第\(i\)项的系数就是\(\begin{Bmatrix}n\\i\end{Bmatrix}\).
于是愉快的将时间变成了\(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\)
类似的，我们可以构造序列\(\{\Delta a_n\}\)的二阶、三阶…\(k\)阶差分序列。 不妨记为\(\{\Delta^2 a_n\},...,\{\Delta^k a_n\}\)
令序列是一个\(p\)次多项式，那么差分表一个很重要且很显然的性质是：$$\forall n>=0, \Delta^{p+1}a_n = 0$$这是由于每次差分都必然会把最高次项消去！
另外一个很重要的性质就是差分表的线性性，即如果\(f_n=k_1g_n+k_2h_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)!}=\tbinom{n}{p}$$
那么根据差分表的线性性，我们就可以得知$$f_n=\sum_{i=0}^p c_i\tbinom{n}{i}$$
由于$$\sum_{k = 0}^n f(k) = \sum_{k = 0}^p c_k {n + 1 \choose k}$$所以利用差分表，我们可以在\(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){n + 1 \choose k}$$
8. 伯努利数
初三再来学一学伯努利数吧。
根据伯努利数的生成函数定义，可知$$\frac{x}{e^{x}-1} = \sum_{i\geq 0} B_{i}*\frac{x^{i}}{i!}$$
由于$$\frac{x}{ex-1}·(ex-1)=x$$
且$$[xn]\frac{x}{ex-1}·(ex-1)=\sum_{i=0}\frac{B_i}{i!}·\frac{1}{(n-i)!}=[n=1]$$
两边都乘上\(n!\)可以得到$$\sum_{i=0}^{n-1}\binom{n}{i}B_i=[n=1]$$
这是伯努利数的一个基本性质。后面会用到。
我们再定义一个多项式\(B_n(t)\)表示$$B_{n}(t) = \sum_{k=0}^{n-1} B_{k} * t{n-k}*C_{n}$$
然后我们发现
注意到后面只有当\(n-i-1=0\)时值为\(1\)，于是$$B_n(t+1)-B_n(t)=n*t^{n-1}$$
然后我们考虑差分，就有$$\sum_{t=0}{n-1}B_k(t+1)-B_k(t)=k·\sum_{i=0}i{k-1}$$即$$B_{k+1}(n+1)=(k+1)·\sum_{i=0}ni^k$$
可得自然数幂和$$\sum_{i=0}nik=\frac{1}{k+1}·\sum_{i=0}kB_in\binom{k+1}{i}$$
问题转化成了求\(B_i\)，注意到它的生成函数定义，事实上我们只需要求$$\sum_{i\ge 0}\frac{x^i}{(i+1)!}$$在模\(x^{k+1}\)的逆元即可。
时间复杂度\(O(KLogK)\)，当然如果递推的话，也是可以轻松做到\(O(k^2)\)的。
9. 拉格朗日插值法
两大作用：
-
快速根据点值逼近原函数.
-
取点值对大于\(n\)时唯一确定\(n\)次多项式.
Example
例如对于$$\sum\limits_{i=1}^ni=\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\)个函数线性结合即可。
Practice
例如我们要求自然数幂和.
各种方法可以证明\(i^k\)的和是\(k+1\)次的, 所以我们只需要给出\(k+2\)个点值表达,就可以求得通项公式.
记\(S(n)=\sum_{i=1}^n i^k\), 则
那么时间复杂度就在预处理\(y_i\)上面了, 利用线筛,可以做到\(O(k)\)级别.
后面的那一部分可以预处理,具体的说,就有:$$S(n)=\sum_{i=1}{k+2}y_ipre[i-1]suf[i+1][(-1)(i-1)!(k+2-i)!]^{-1}$$
\(Code\)
时间复杂度：\(O(\frac{k}{ln_k}log_2k)=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);
}
