求解自然数幂和的若干种方法

Posted 2021-02-15 pro-king

问题的引入

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

1. 循环

四年级应该会循环了。

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

2. 快速幂

五年级学了快速幂之后就能做到(O(nlog_2k))

请不要小看这个算法。有时候在特定的情况下（例如(n)很小，或(1 ightarrow 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}^nsum_{j=0}^{k-1}C_k^ji^j\ =sum_{i=0}^{k-1}C_k^isum_{j=1}^n j^i\ =sum_{i=0}^{k-1}C_k^i S(i) ][ herefore (n+1)^{k+1}-1=sum_{i=0}^kC_{k+1}^iS(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)=f(frac{n}{2},k)+sum_{j=0}^kC_k^j*f(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)]

但事实上，我们求解自然数幂和需要用到的是它的原始定义：
[x^{ndownarrow}=xcdot (x-1)cdot (x-2)cdots (x-n+1)=sum_{k=0}^ns_s(n,k)cdot x^k][x^{nuparrow}=xcdot (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^{nuparrow}=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^{mdownarrow}-sum_{k=0}^{m-1}S_s(m,k)·n^k]

我们继续推，发现下降幂的和是可以写成一个组合数的形式的，比方说[sum_{i=m}^ni^{mdownarrow}=sum_{i=m}^nfrac{i!m!}{(i-m)!m!}=m!sum_{i=m}^ninom{i}{m}=m!inom{n+1}{m+1}]

而后面那一坨式子也是可以化简的，比方说[sum_{i=0}^nsum_{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)个集合，且集合非空的方案数。

基本性质是[egin{Bmatrix}n\mend{Bmatrix}=egin{Bmatrix}n-1\m-1end{Bmatrix}+mcdot egin{Bmatrix}n-1\mend{Bmatrix}]

考虑它的通项公式，可以先把所有集合标号，最后除以集合的阶乘即可，那么考虑容斥，枚举非空集合个数(i)，可以得到[egin{Bmatrix}n\mend{Bmatrix}=frac 1 {m!}sum_{i=0}^m(-1)^iinom mi(m-i)^n]

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

显然！！[i^k=sum_{j=0}^kleft{egin{array}{c}{k}\{j}end{array} ight}i^{jdownarrow}]

继续推导[sum_{i=1}^ni^k=sum_{i=1}^nsum_{j=0}^kleft{egin{array}{c}{k}\{j}end{array} ight}i^{jdownarrow}][=sum_{i=1}^nsum_{j=0}^kleft{egin{array}{c}{k}\{j}end{array} ight} j!left(egin{array}{c}{i}\{j}end{array} ight)][=sum_{j=0}^kleft{egin{array}{c}{k}\{j}end{array} ight} j!sum_{i=j}^nleft(egin{array}{c}{i}\{j}end{array} ight)][=sum_{j=0}^kleft{egin{array}{c}{k}\{j}end{array} ight} j!inom{n+1}{j+1}]

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

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

然而事实上，我们来看看斯特林数的通项公式：[egin{Bmatrix}n\mend{Bmatrix}=frac 1 {m!}sum_{i=0}^m(-1)^iinom mi(m-i)^n]

一拼凑，咦~[egin{Bmatrix}n\mend{Bmatrix}=sum_{i=0}^mfrac{(-1)^i}{i!}cdotfrac{(m-i)^n}{(m-i)!}]

这原来可以写成形如(sum_{i=0}^mf(i)*g(m-i))的卷积形式。于是第一个多项式的第(i)项系数是(frac {(-1)^i}{i!})，另一个多项式的第(i)项系数是(frac {i^n}{i!})，卷积后第(i)项的系数就是(egin{Bmatrix}n\iend{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_1Delta^p g_n + k_2Delta^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)!}= binom{n}{p}]

那么根据差分表的线性性，我们就可以得知[f_n=sum_{i=0}^p c_i binom{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_{igeq 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}inom{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} group 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}inom{k+1}{i}]

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

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

9. 拉格朗日插值法

两大作用：

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

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

Example

例如对于[sumlimits_{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{prodlimits_{j eq i}(x-x_j)}{prodlimits_{j eq 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_iprod_{j=1,i eq j}^{k+2}frac {n-x_j}{x_i-x_j}=sum_{i=1}^{k+2}y_ifrac {prod_{j=1,i eq j}^{k+2}(n-j)}{prod_{j=1,i eq 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
", (Solve(r) - Solve(l - 1) + Mo) % Mo);
}