自然数幂和

自然数幂和问题

给定 \(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\) 次的多项式。

代码：

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