题意:
给定n(n<=1e9)和k(k<=1e6),求出\(\sum_{i=1}^n i^k\)
拉格朗日的经典题目
首先题目中提到,\(\sum_{i=1}^n i^k\)可以转化为以\(n\)为未知数的\(k+1\)次多项式\(y=a_1x^{k+1}+a_2x^k+a_3x^{k-1}+...+a_{k+2}x^0\)。
那么实际上就变成了,当x等于n时求这个多项式的值
根据拉格朗日差值法,需要k+2个点才能确定多项式,但是时间是\(O(k^2)\),完不成。
但是,由于给定的点的x可以为连续的值,所以可以优化为\(O(n)\)
方法:
\[y=\sum_{i=1}^{k+2}y_i\prod_{j=1,j\neq i}^{k+2}(x-x_j)/(x_i-x_j)
\]
我们可以看到
\(y_i\)可以通过\(y_{i-1}+i^k\)递推得到复杂度为\(O(\log k)\)
用前缀和的思想可以处理出pre[i]为\(\prod_{j=1}^i(x-x_j)\)
同样用相同的思想处理处suc[i]为\(\prod_{j=i}^{k+2}(x-x_j)\)
这样\(\prod_{j=1,j\neq i}^{k+2}(x-x_j)\)就可以用pre[i-1]*suc[i+1]来\(O(1)\)
而各个点的x为连续的,所以\(\prod_{j=1,j\neq i}^{k+2}(x_i-x_j)\)中j<i的部分,乘积刚好是i-j的阶乘。
同样j>i的部分也是k+2-i的阶乘,但是这一部分的正负需要判断!
所以提前处理阶乘可以实现\(O(1)\)求出。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=1e9+7;
const int maxk=1e6+10;
int n,k,ans;
int y[maxk];
int jc[maxk],pre[maxk],suc[maxk];
int qpow(int a, int b) {
int ans = 1;
for(; b >= 1; b >>= 1, a = 1ll * a * a % mod)
if(b & 1) ans = 1ll * ans * a % mod;
return ans;
}
int main()
{
scanf("%d%d",&n,&k);
jc[0]=1;
for(int i=1;i<=k+2;++i)jc[i]=1ll*jc[i-1]*i%mod;
pre[0]=1;
for(int i=1;i<=k+2;++i)pre[i]=1ll*pre[i-1]*(n-i)%mod;
suc[k+3]=1;
for(int i=k+2;i>0;--i)suc[i]=1ll*suc[i+1]*(n-i)%mod;
y[1]=1;
for(int i=2;i<=k+2;++i)y[i]=(y[i-1]+qpow(i,k))%mod;
for(int i=1;i<=k+2;++i)
{
int fz=1ll*pre[i-1]*suc[i+1]%mod*y[i]%mod;
int fm=1ll*jc[i-1]*jc[k+2-i]*(((k+2-i)&1)?-1:1)%mod;
ans=(ans+1ll*fz*qpow(fm,mod-2)%mod)%mod;
ans=(1ll*ans+mod)%mod;
}
cout<<ans<<endl;
return 0;
}