之前写过 自然数方幂和公式， 这次写的目的是因为上次是从数学上完美的解决了这个问题，这次我们要从计算上完美的解决这个问题，当然这归功于我看到的一份 sgtlaugh 的代码。经过解读体会到其中的奥秘，特此记录。一句话，简直不敢相信。

如果有人说他能在 $O(k)$ 时空复杂度求解 $\sum_{i=1}^n i^k$，你肯定会说这怎么可能别忽悠我了，那我只能说，因为你没看过这篇博文。

首先预处理，$O(k)$ 复杂度 求 $\sum_{i=1}^n i^k$ 其中 $n \leq k$

我之前一直以为要用 $k \log k$ 的复杂度才能解决这个问题，其实我们只需对所有素数 $p$ 计算 $p^k$ 即可。对于一般的 $i$ 我们先预处理 其最小素因子 $sp[i]$。计算 $sp[i]^k \cdot (i/sp[i])^k$ 即可（具体可见最后代码）。由于素数的阶为 $O(\frac{k}{\log k})$ 因此整个复杂度即为 $O(k)$。

再由 Lagrange 插值多项式得出最终答案

因为我们知道 $\sum_{i=1} ^n i^k$ 一定是一个关于 $n$ 的次数为 $k+1$ 的多项式。因此，我们只需计算其在 $0,\cdots,k+1$ 上的取值，用 Lagrange 插值多项式即可知道答案。

对于一个次数不超过 $n$ 的多项式 $f(x)$，其在不同位置 $x_0,\cdots,x_n$ 的取值唯一决定了这个多项式:

具体到本问题，我们取 $x=n, m=k+1, x_i=i$ 那么

例题：Codeforces 622F

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#include <bits/stdc++.h>
typedef long long LL;
using namespace std;
const int N = 1e6+6;
const LL mod = 1e9+7;
int sp[N],p[N];
int spf(){ // samllest prime factor
    int cnt=0;p[cnt++]=2;
    for(int i=2;i<N;i+=2)   sp[i]=2;
    for(int i=1;i<N;i+=2)   sp[i]=i;
    for(int i=3;i<N;i+=2){
        if(sp[i]==i)    p[cnt++] = i;
        for(int j = 1;j < cnt && p[j]<=sp[i] && i * p[j] < N; ++j) {
            sp[i * p[j]] = p[j];
        }
    }
    return cnt;
}
LL pow_mod(LL x,LL n,LL p){
    LL r=1;
    while(n){
        if(n&1) r=r*x%p;
        n>>=1;  x=x*x%p;
    }
    return r;
}
LL inv[N],AP[N],AS[N],f[N];
LL getpowsum(LL n,int k){ // mod > k
    if(k==0)    return n%mod;
    if(p[0]!=2) spf();
    int nk=k+1;
    f[0]=0;f[1]=1;
    for(int i=2;i<=nk;++i){
        if(sp[i]==i)    f[i]=pow_mod(i,k,mod);
        else    f[i]=f[sp[i]]*f[i/sp[i]]%mod;
    }
    for(int i=1;i<=nk;++i){
        f[i]+=f[i-1];
        if(f[i]>=mod)   f[i]-=mod;
    }
    if(n<=nk)    return f[n];
    LL tmp = 1;
    for(int i=2;i<=nk;++i)    tmp=tmp*i%mod;
    inv[nk] = pow_mod(tmp,mod-2,mod);
    for(int i=nk-1;i>=0;--i) inv[i]=inv[i+1]*(i+1)%mod;
    AP[0]=AS[nk]=1;
    for(int i=1;i<=nk;++i)   AP[i]=AP[i-1]*(n+1-i)%mod;
    for(int i=nk-1;i>=0;--i) AS[i]=AS[i+1]*(n-i-1)%mod;
    LL res = 0;
    for(int i=1;i<=nk;++i){ // because f(i)=0
        LL x = f[i]*AP[i]%mod*AS[i]%mod*inv[i]%mod*inv[nk-i]%mod;
        if((nk-i)&1) res-=x; // be careful
        else        res+=x;
    }
    return (res%mod+mod)%mod;
}
int main(){
    LL n;
    int k;
    while(cin>>n>>k){
        cout<<getpowsum(n,k)<<endl;
    }
    return 0;
}

实际上我们可以不求 $\mod p$ 后的答案，利用大数类得到标准答案，但是这时因为数字实在太大，每次乘法的用时过大，因此仅适合 $k<n$ 的情况

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
cpp_int f[N];
cpp_int getpowsum(LL n,int k){
    if(k==0)    return cpp_int(n);
    if(p[0]!=2) spf();
    int nk=2*k+1;
    f[0]=0;f[1]=1;
    for(int i=2;i<=nk+1;++i){
        if(sp[i]==i)    f[i]=pow(cpp_int(i),k);
        else    f[i]=f[sp[i]]*f[i/sp[i]];
    }
    for(int i=1;i<=nk;++i)  f[i]+=f[i-1];
    if(n<=nk)   return f[n];
    cpp_int res = 0,tl=1,tr=1;
    for(int i=nk-1;i>=0;--i)    tr=tr*(n-i-1)/(nk-i);
    for(int i=0;i<=nk;++i){
        if((nk-i)&1)    res -= f[i]*tl*tr;
        else            res += f[i]*tl*tr;
        tl = tl*(n-i)/(i+1);
        tr = tr*(nk-i)/(n-i-1);
    }
    return res;
}

其实如果我们知道最终的上界，求出多个 $\mod p$ 后的答案，再用中国剩余定理貌似很不错。

该方法可以推广成求 $\sum_{i=1}^n f(i)^k$，其中 $f(x)$ 是多项式。具体分析即可

这种情况一般很难再做到 $O(k)$ 时间复杂度，而变成了 $O(k \log k) \deg f$ 复杂度。

如果求所有 $f(k) = \sum_{i = 1}^n i^k, k = 0, \cdots, m$，可以在 $O(m \log m)$ 复杂度求出

$\frac{f(k)}{k!}$ 是 $\sum_{i = 1}^n e^{ix} = e^x \frac{e^{nx} - 1}{e^x - 1}$ 的 $x^k$ 的系数。注意到 $e^x - 1$ 的常数项系数为 0，所以不可逆，要分子分母同时除以 $x$，所以求它们的时候要多预算一位。

这里用幂生成函数而非常规生成函数是因为 $\sum_{i = 1}^n \frac{1}{1 - i x}$ 不好计算，维护分子分母分治可以在 $O(n \log^2 n)$ 计算出结果。