自然数幂和

自然数幂和问题

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

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

拉格朗日插值解法

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

$$\begin{array}{r}f(x)=\sum _{i=1}^n{y}_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}\end{array}$$

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

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

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

$$\begin{array}{rl}\sum _{i=1}^n{i}^k& =\sum _{i=1}^n\sum _{j=0}^{k}S(k,j)i^{\underset{―}{j}}\\ & =\sum _{j=0}^{k}S(k,j)\sum _{i=1}^n{i}^{\underset{―}{j}}\end{array}$$

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

$$\begin{array}{r}\mathrm{\Delta }{i}^{\underset{―}{j}}=j\times i^{\underset{―}{j-1}}\end{array}$$

则原式可以化为：

$$\begin{array}{r}\sum _{j=0}^{k}S(k,j)\sum _{i=1}^n{i}^{\underset{―}{j}}=\sum _{j=0}^{k}S(k,j)\sum _{i=1}^n\frac{\mathrm{\Delta }{i}^{\underset{―}{j+1}}}{j}\end{array}$$

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

代码：

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;
}