自然数幂和问题
给定 \(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;
}