自然数幂和问题

给定 $n,k$，求 $\sum_{i=1}^n i^k$ 的值。

$n\le 10^9,k\le 10^6$

拉格朗日插值解法

给出拉格朗日插值的式子：

$$
\begin{aligned}
f(x)=\sum\limits_{i=1}^n y_i\prod\limits_{j\not =i}\frac{x-x_j}{x_i-x_j}
\end{aligned}
$$

容易看出计算 $f(x)$ 的时间复杂度为 $O(n^2)$，但是如果 $x_i$ 连续，且 $x$ 为整数，那右边分式分母可以写成阶乘的形式。可以在 $O(n)$ 的时间复杂度内求解。

对于原问题，显然 $n$ 应该是变量，现在证明原问题是关于 $n$ 的 $k+1$ 次多项式。

根据第二类斯特林数的经典结论，我们有：

$$
\begin{aligned}
\sum\limits_{i=1}^ni^k&=\sum\limits_{i=1}^n\sum\limits_{j=0}^kS(k,j)i^{\underline{j}}\
&=\sum\limits_{j=0}^k S(k,j) \sum\limits_{i=1}^n i^{\underline{j}}
\end{aligned}
$$

设 $\Delta f(x)=f(x+1)-f(x)$，即给出有限微积分的形式，则有：

$$
\begin{aligned}
\Delta i^{\underline{j}}=j\times i^{\underline{j-1}}
\end{aligned}
$$

则原式可以化为：

$$
\begin{aligned}
\sum\limits_{j=0}^k S(k,j) \sum\limits_{i=1}^n i^{\underline{j}}=\sum\limits_{j=0}^k S(k,j)\sum\limits_{i=1}^n\frac{\Delta i^{\underline{j+1}}}{j}
\end{aligned}
$$

其中，后者是一个 $k+1$ 次的多项式。

代码：

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
int x[N],y[N],n,k,sum;
int jie[N],invjie[N];

inline int ksm(int a,int b,int mod){
    int res=1;while(b){if(b&1)res=1ll*res*a%mod;a=1ll*a*a%mod;b>>=1;}return res;
}
inline int inv(int x){return ksm(x,mod-2,mod);}
inline int Leg(int h){
    sum=1;rep(i,1,k+2) sum=1ll*sum*(h-x[i])%mod;
    rep(i,1,k+2) x[i]=sum*inv(h-x[i])%mod;
    int ans=0;
    rep(i,1,k+2){
        ans=(ans+y[i]*x[i]%mod*invjie[i-1]%mod*invjie[k+2-i]%mod*(((k+2-i)&1)?(mod-1):1)%mod)%mod;
    }
    return ans;
}

signed main(){
    read(n);read(k);jie[0]=1;
    rep(i,1,k+2) jie[i]=jie[i-1]*i%mod;invjie[k+2]=inv(jie[k+2]);
    dec(i,0,k+1) invjie[i]=invjie[i+1]*(i+1)%mod;
    rep(i,1,k+2){
        x[i]=i;y[i]=(ksm(i,k,mod)+y[i-1])%mod;
    }
    if(n<=k+2){
        printf("%lld\n",y[n]);
    }
    else{
        int ans=Leg(n);
        printf("%lld\n",ans%mod);
    }
    return 0;
}