問題的引入

給定\(n,k\)求\[\sum_{i=1}^ni^k\]

1. 循環

四年級應該會循環了。

能做到\(O(nk)\)的優秀時間復雜度。

2. 快速冪

五年級學了快速冪之後就能做到\(O(nlog_2k)\)

請不要小看這個算法。有時候在特定的情況下（例如\(n\)很小，或\(1\rightarrow n\)的距離變得很小時），這個復雜度真的很優秀。

3. 差分法

六年級應該知道差分和二項式定理了。那麽：\[(a+1)^k-a^k=\sum_{i=0}^{k-1}C_k^ia^i\]

於是：
\[ (n+1)^k-1 =\sum_{i=1}^n (i+1)^k-i^k \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =\sum_{i=1}^n\sum_{j=0}^{k-1}C_k^ji^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) \]

把\(i=k\)時移項，可以得到\[(k+1)S(k)=(n+1)^{k+1}-1-\sum_{i=0}^{k-1}C_{k+1}^iS(i)\]

所以\[S(k)=\frac{(n+1)^{k+1}-1-\sum_{i=0}^{k-1}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}\)我們計算的時間復雜度都是\(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)\]

但事實上，我們求解自然數冪和需要用到的是它的原始定義：
\[x^{n\downarrow}=x\cdot (x-1)\cdot (x-2)\cdots (x-n+1)=\sum_{k=0}^ns_s(n,k)\cdot x^k\]\[x^{n\uparrow}=x\cdot (x+1)\cdot (x+2)\cdots(x+n-1)=\sum_{k=0}^ns_u(n,k)\cdot x^k\]

這裏需要註意，第一類斯特林數根據定義分成了有符號\(S_s\)和無符號\(S_u\)兩種。事實上，我們可以很輕松的從這個原始定義推出它的組合意義。

因為\[\sum_{k=0}^nS_u(n,k)·x^k=x^{n\uparrow}=x^{(n-1)\uparrow}·(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}^nS_s(n,k)x^k\]

於是我們可以得到一個顯然的式子是：\[n^m=n^{m\downarrow}-\sum_{k=0}^{m-1}S_s(m,k)·n^k\]

我們繼續推，發現下降冪的和是可以寫成一個組合數的形式的，比方說\[\sum_{i=m}^ni^{m\downarrow}=\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}^{m-1}S_s(m,k)·i^k=\sum_{k=0}^{m-1}S_s(m,k)\sum_{i=0}^ni^k\]

發現後面那條式子\(\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\]

接下來繼續推導自然數冪和。

顯然！！\[i^k=\sum_{j=0}^k\left\{\begin{array}{c}{k}\\{j}\end{array}\right\}i^{j\downarrow}\]

繼續推導\[\sum_{i=1}^ni^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!\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}{e^x-1}·(e^x-1)=x\]

且\[[x^n]\frac{x}{e^x-1}·(e^x-1)=\sum_{i=0}^{n-1}\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}^{k}\]

然後我們發現
\[B_n(t+1)-B_n(t)=\sum_{k=0}^{n-1} B_{k}*\lgroup(t+1)^{n-k} - t^{n-k}\rgroup C_{n}^{k}\]\[=\sum_{k=0}^{n-1} B_{k}*(\sum_{i=0}^{n-k-1} C_{n-k}^{i} * t^{i})* C_{n}^{k}\]\[=\sum_{k=0}^{n-1} B_{k}*(\sum_{i=0}^{n-k-1} C_{n-k}^{i} * t^{i}*C_{n}^{k})\]\[=\sum_{i=0}^{n-1}B_k*\sum_{i=0}^{n-k-1}\frac{n!t^i}{i!k!(n-k-i)!}\]\[=\sum_{k=0}^{n-1} B_{k}*(\sum_{i=0}^{n-k-1} C_{n-i}^{k} * t^{i}*C_{n}^{i})\]\[=\sum_{i=0}^{n-1} C_{n}^{i}*t^{i}*\sum_{k=0}^{n-1-i} B_{k}*C_{n-i}^{k}\]

註意到後面只有當\(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}^{n-1}i^{k-1}\]即\[B_{k+1}(n+1)=(k+1)·\sum_{i=0}^ni^k\]

可得自然數冪和\[\sum_{i=0}^ni^k=\frac{1}{k+1}·\sum_{i=0}^kB_in^{k+1-i}\binom{k+1}{i}\]

問題轉化成了求\(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\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\)個函數線性結合即可。

\[f_i(x)=\frac{\prod\limits_{j\neq i}(x-x_j)}{\prod\limits_{j\neq i}(x_i-x_j)}*y_i\]\[g(x)=\sum_{i=0}^nf_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\neq j}^{k+2}\frac {n-x_j}{x_i-x_j}=\sum_{i=1}^{k+2}y_i\frac {\prod_{j=1,i\neq j}^{k+2}(n-j)}{\prod_{j=1,i\neq j}^{k+2}(i-j)}\]

那麽時間復雜度就在預處理\(y_i\)上面了, 利用線篩,可以做到\(O(k)\)級別.

後面的那一部分可以預處理,具體的說,就有:\[S(n)=\sum_{i=1}^{k+2}y_ipre[i-1]suf[i+1][(-1)^{k+2-i}(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);
}

求解自然數冪和的若幹種方法