浅谈高次幂求和

· · 算法·理论

须知

参考了 oi-wiki 的 this（里面也有许多本帖未提及的证明）插值部分参考了 this。

这是本帖用到的所有不常见的数学符号：

S_k(n)=\sum_{i=1}^{n}i^k=1^k+2^k+3^k+\dots+n^k \dbinom{n}{m}=\frac{n!}{m!(n-m)!}

引入

我们都知道 0 到 3 次的求和公式：

S_0(n)=n S_1(n)=\frac{n(n+1)}{2} S_2(n)=\frac{n(n+1)(2n+1)}{6} S_3(n)=\frac{n^2(n+1)^2}{4}

即：

S_0(n)=n S_1(n)=\frac{1}{2}n^2+\frac{1}{2}n S_2(n)=\frac{1}{3}n^3+\frac{1}{2}n^2+\frac{1}{6}n S_3(n)=\frac{1}{4}n^4+\frac{1}{2}n^3+\frac{1}{4}n^2

那么 k 次的求和公式呢？

其实可以用伯努利数计算得到。

假设我们已知零，一次的求和公式，现在要求二次求和公式。

(x-1)^3=x^3-3x^2+3x-1

移项：

-x^3+(x-1)^3=-3x^2+3x-1 x^3-(x-1)^3=3x^2-3x+1

将 x= 不同数代入：

1^3-0^3=3\times 1^2-3\times 1+1 2^3-1^3=3\times 2^2-3\times 2+1 3^3-2^3=3\times 3^2-3\times 3+1 \dots n^3-(n-1)^3=3n^2-3n+1

全部相加：

n^3=3\sum_{i=1}^{n}i^2-3\sum_{i=1}^{n}i+n

即：

n^3=3S_2(n)-3S_1(n)+S_0(n)

移项：

3S_2(n)=3S_1(n)-S_0(n)+n^3

我们已知 S_0(n)，S_1(n)，代入。

3S_2(n)=3\frac{n^2+n}{2}-n+n^3 3S_2(n)=\frac{2n^3+3n^2+n}{2} S_2(n)=\frac{2n^3+3n^2+n}{6}

即：

S_2(n)=\frac{1}{3}n^3+\frac{1}{2}n^2+\frac{1}{6}n

对于多次方和，可以使用相同的方法得到。

现在将它们排列在一起：

S_0(n)=\frac{1}{1}n S_1(n)=\frac{1}{2}n^2+\frac{1}{2}n S_2(n)=\frac{1}{3}n^3+\frac{1}{2}n^2+\frac{1}{6}n S_3(n)=\frac{1}{4}n^4+\frac{1}{2}n^3+\frac{1}{4}n^2+0\times n S_4(n)=\frac{1}{5}n^5+\frac{1}{2}n^4+\frac{1}{3}n^3+0\times n^2-\frac{1}{30}n S_5(n)=\frac{1}{6}n^6+\frac{1}{2}n^5+\frac{5}{12}n^4+0\times n^3-\frac{1}{12}n^2+0\times n

可以将系数取出来作为一个矩阵：

\begin{bmatrix} \frac{1}{1}&&&&&\\ \frac{1}{2}&\frac{1}{2}&&&&\\ \frac{1}{3}&\frac{1}{2}&\frac{1}{6}&&&\\ \frac{1}{4}&\frac{1}{2}&\frac{1}{4}&0&&\\ \frac{1}{5}&\frac{1}{2}&\frac{1}{3}&0&\frac{1}{30}&\\ \frac{1}{6}&\frac{1}{2}&\frac{5}{12}&0&\frac{1}{12}&0\\ \end{bmatrix}

注意到（注意力惊人）：

X_{i,j}=X_{0,j}\times X_{i,0}\times \dbinom{i+1}{j}

也就是说，第 i 行第 j 列的系数只与它这一行最左边的数和这一列最上面的数有关。

最左边的数的我们可以发现，它为 \dfrac{1}{i+1}。

所以现在要求 k 次方和公式，只需知道每列最上面（即每行最左边）的数就可以了。

事实上，它就是伯努利数。

伯努利数

数列 B 表示伯努利数，其前几项为：

1,\frac{1}{2},0,\frac{1}{6},0,-\frac{1}{30},0,\frac{1}{42}\cdots

可以发现 B_{2k+1}=0(k\in\mathbb{N}^+)。

那它怎么求呢？

这里不加证明的给个公式：

\sum_{k=0}^{n}\dbinom{n}{k}B_k=nB_n(n\ge 1)

再根据上面的推导，我们便可以求出 k 次方和公式了：

\sum_{i=1}^{n}i^k=\frac{1}{k+1}\sum_{i=0}^k\dbinom{k+1}{i}B_i(n+1)^{k-i+1}

代码实现：

由于答案太大，通常加上取模，用逆元代替除法：

#include<bits/stdc++.h>
#define inv(x) power(x,mod-2)
using namespace std;
constexpr int mod=998244353;
int n;
vector<int>ans;
int vis[5005][5005];
int power(int x,int y){
    int res=1;
    while(y){
        if(y&1) res=1ll*res*x%mod;
        x=1ll*x*x%mod,y>>=1;
    }
    return res;
}
int C(int n,int m){
    if(vis[n][m]) return vis[n][m];
    return vis[n][m]=!m||n==m?1:(C(n-1,m-1)+C(n-1,m))%mod;
}
vector<int>get(int n){
    vector<int>ans(n+1,0);
    ans[0]=1;
    for(int i=1;i<=n;i++){
        if(i>1&&(i&1)) continue;
        for(int j=0;j<i;j++) ans[i]=(1ll*ans[j]*C(i+1,j)+ans[i])%mod;
        ans[i]=1ll*ans[i]*inv(i+1)%mod,ans[i]=-ans[i];
    }
    return ans;
}
int main(){
    scanf("%d",&n),ans=get(n);
    for(auto p:ans) printf("%d\n",p);
    return 0;
}

求出伯努利数，我们便可以求出 k 次方和了：

#include<bits/stdc++.h>
#define inv(x) power(x,mod-2)
using namespace std;
constexpr int mod=998244353;
int k,n;
vector<int>ans;
int vis[5005][5005];
int power(int x,int y){
    int res=1;
    while(y){
        if(y&1) res=1ll*res*x%mod;
        x=1ll*x*x%mod,y>>=1;
    }
    return res;
}
int C(int n,int m){
    if(vis[n][m]) return vis[n][m];
    return vis[n][m]=!m||n==m?1:(C(n-1,m-1)+C(n-1,m))%mod;
}
vector<int>get(int n){
    vector<int>ans(n+1,0);
    ans[0]=1;
    for(int i=1;i<=n;i++){
        if(i>1&&(i&1)) continue;
        for(int j=0;j<i;j++) ans[i]=(1ll*ans[j]*C(i+1,j)+ans[i])%mod;
        ans[i]=1ll*ans[i]*inv(i+1)%mod,ans[i]=-ans[i];
    }
    return ans;
}
int sum(int k,int n){
    vector<int>B=get(k);int ans=0;
    for(int i=0;i<=k;i++) ans=((__int128)C(k+1,i)*B[i]*power(n+1,k-i+1)*inv(k+1)+ans)%mod;
    return ans;
}
int main(){
    scanf("%d%d",&k,&n);
    printf("%d",sum(k,n));
    return 0;
}

值得一提的是，这个做法的时间瓶颈在于求伯努利数，时间复杂度为 \Theta(k^2)，但我们使用插值可以做到 \Theta(k) 的时间复杂度。

插值做法

我们都知道插值公式：

f(x)=\sum_{i=1}^{n}y_{i}\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}

但直接插值显然是 \Theta(k^2) 的，因此我们要想办法优化。

在这里，横坐标显然都是连续整数，因此可以优化到 \Theta(k)。

因为横坐标是连续整数，所以可以假设 x_i=i，那么上面的公式变成：

f(x)=\sum_{i=1}^{n}y_{i}\prod_{j\ne i}\frac{x-j}{i-j}

先考虑分子：

\prod_{j\ne i}(x-j)=\frac{\prod\limits_{j=1}^{n+1}(x-j)}{x-i}

这就可以预处理 x-i 的后缀积，做到 \Theta(1) 求分子。

而分母就将它拆成阶乘，使其不含有 j：

\prod_{j\ne i}(i-j)=(-1)^{n-i+1}(i-1)!(n-i+1)!

于是我们就可以通过预处理阶乘和逆元来做到 \Theta(1) 求分母。

全部代入就可以得到 \Theta(n) 的插值公式:

f(x)=\sum_{i=1}^{n+1}(-1)^{n-i+1}y_{i}\frac{\prod\limits_{j=1}^{n+1}(x-j)}{(i-1)!(n-i+1)!(x-i)}

现在我们回到问题。显然，S_k(n) 是一个 k+1 次多项式，因此我们可以筛出 1^i,\cdots,(k+2)^i 的值然后进行插值。

其实这个问题还可以组合数学相关知识直接推出答案式子，感兴趣的可以去最上面须知的 oi-wiki 链接查看。

代码实现：

#include<bits/stdc++.h>
using namespace std;
constexpr int mod=998244353;
int n,k,cnt,ans,flag[1000005],prime[1000005],f[1000005];
int pre[1000005],suf[1000005],fac[1000005],inv[1000005];
int power(int x,int y){
    int res=1;
    while(y){
        if(y&1) res=1ll*res*x%mod;
        x=1ll*x*x%mod,y>>=1;
    }
    return res;
}
void sieve(int x){
    f[1]=1;
    for(int i=2;i<=x;i++){
        if(!flag[i]) prime[++cnt]=i,f[i]=power(i,k);
        for(int j=1;j<=cnt;j++){
            if(1ll*i*prime[j]>x) break;
            flag[i*prime[j]]=true;
            f[i*prime[j]] =1ll*f[i]*f[prime[j]]%mod;
            if(!(i%prime[j])) break;
        }
    }
    for(int i=2;i<=x;i++) f[i]=(f[i]+f[i-1])%mod;
}
int main() {
    scanf("%d%d",&k,&n),sieve(k+2);
    if(n<=k+2) return printf("%d",f[n]),0;
    pre[0]=suf[k+3]=fac[0]=inv[0]=fac[1]=inv[1]=1;
    for(int i=1;i<=k+2;i++) pre[i]=1ll*pre[i-1]*(n-i)%mod;
    for(int i=k+2;i>=1;i--) suf[i]=1ll*suf[i+1]*(n-i)%mod;
    for(int i=2;i<=k+2;i++){
        fac[i]=1ll*fac[i-1]*i%mod;
        inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    }
    for(int i=2;i<=k+2;i++) inv[i]=1ll*inv[i]*inv[i-1]%mod;
    for(int i=1;i<=k+2;i++){
        int x=1ll*pre[i-1]*suf[i+1]%mod;
        int y=1ll*inv[i-1]*inv[k-i+2]%mod;
        int neg=((k+i)&1)?-1:1;
        ans=(ans+1ll*(x*neg+mod)%mod*y%mod*f[i]%mod)%mod;
    }
    printf("%d",ans);
    return 0;
}