51nod 1258 序列求和 V4 拉格朗日插值法求自然数幂和

题意

T(n) = n^k，S(n) = T(1) + T(2) + …… T(n)。给出n和k，求S(n)。
例如k = 2，n = 5，S(n) = 1^2 + 2^2 + 3^2 + 4^2 + 5^2 = 55。
由于结果很大，输出S(n) Mod 1000000007的结果即可。
1 <= T <= 500，1 <= N <= 10^18, 1 <= K <= 50000

分析

拉格朗日插值法：
如果知道一个n次多项式f(x)上的n+1个点(x0,y0),(x1,y1),…,(xn,yn)，现在我们知道x’，要求y’。
那么我们可以构造一个多项式$f(x)=b_0(x)*y_0+b_1(x)*y_1+...+b_n(x)*y_n$其中$b_k(x)$均为n次多项式。
要使这个多项式满足这n+1个点，则有$b_k(x_k)=1,b_k(x_{k'})=0(k!=k')$，那么$x_{k'}$必然是$b_k(x)$的一个根。
不难发现有

$$b_k(x)=\frac{\prod_{i=0,i!=k}^n(x-x_i)}{\prod_{i=0,i!=k}^n(x_k-x_i)}$$
那么我们就可以通过把x’代进这条式子里面来求出y’。

再回到这题，设$ans(x)=\sum_{i=1}^xi^k$，那么显然ans必然是一个k+1次的多项式。我们可以先通过快速幂把$ans(i)(1<=i<=k+2)$求出，然后通过拉格朗日乘数法来求出$ans(n)$。注意到每一项分子和分母都是两段连续的数相乘，那么我们就可以在较短的时间内把该式子算出来。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;

typedef long long LL;

const int N=50005;
const int MOD=1000000007;

int k,a[N],ny[N],pow[N],p[N],q[N];
LL n;

int ksm(int x,int y)
{
    int ans=1;
    while (y)
    {
        if (y&1) ans=(LL)ans*x%MOD;
        x=(LL)x*x%MOD;y>>=1;
    }
    return ans;
}

int solve(LL x)
{
    if (x<=k+2) return a[x];
    x%=MOD;
    int w=ksm(pow[k+2],MOD-2);LL ans=0;
    p[0]=q[k+3]=1;
    for (int i=1;i<=k+2;i++) p[i]=(LL)p[i-1]*(x-i)%MOD;
    for (int i=k+2;i>=1;i--) q[i]=(LL)q[i+1]*(x-i)%MOD;
    for (int i=1;i<=k+2;i++)
    {
        (ans+=(LL)a[i]*w%MOD*p[i-1]%MOD*q[i+1])%=MOD;
        w=(LL)w*(i-k-2)%MOD*ny[i]%MOD;
    }
    ans+=ans<0?MOD:0;
    return ans;
}

int main()
{
    int T;
    scanf("%d",&T);
    for (int i=1;i<=50002;i++) ny[i]=ksm(i,MOD-2);
    pow[1]=1;
    for (int i=2;i<=50002;i++) pow[i]=(LL)pow[i-1]*(1-i)%MOD;
    while (T--)
    {
        scanf("%lld%d",&n,&k);
        for (int i=1;i<=k+2;i++) a[i]=a[i-1]+ksm(i,k),a[i]-=a[i]>=MOD?MOD:0;
        printf("%d\n",solve(n));
    }
    return 0;
}