cf F. The Sum of the k-th Powers (拉格朗日插值法求自然数幂和）模板

题目链接：哆啦A梦传送门

题意：求自然数幂和 ：   (n<1e9，k<1e6)。

题解：直接拉格朗日插值法就求出来了。

设  ，我们要想用拉格朗日插值法，就首先要知道f(n)函数的最高次幂才行。这里我们用差分就很显然了。

f(n)-f(n-1)=n^k，很显然f(n)的最高次幂为k+1，为什么呢？ f '(n)=f(n)-f(n-1)=n^k，显然原函数最高次幂为k+1。

 

拉格朗日插值法参考博客：

https://blog.csdn.net/lvzelong2014/article/details/79159346

维基

https://www.cnblogs.com/ECJTUACM-873284962/p/6833391.html

 

杜教模板：

/// 注意mod，使用前须调用一次 polysum::init(int M);
namespace polysum {
    #define rep(i,a,n) for (int i=a;i<n;i++)
    #define per(i,a,n) for (int i=n-1;i>=a;i--)
    const int D=1010000;///可能需要用到的最高次
    LL a[D],f[D],g[D],p[D],p1[D],p2[D],b[D],h[D][2],C[D];
    LL powmod(LL a,LL b){LL res=1;a%=mod;assert(b>=0);for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;}

     ///函数用途：给出数列的（d+1）项，其中d为最高次方项
    ///求出数列的第n项，数组下标从0开始
    LL calcn(int d,LL *a,LL n) { /// a[0].. a[d]  a[n]
        if (n<=d) return a[n];
        p1[0]=p2[0]=1;
        rep(i,0,d+1) {
            LL t=(n-i+mod)%mod;
            p1[i+1]=p1[i]*t%mod;
        }
        rep(i,0,d+1) {
            LL t=(n-d+i+mod)%mod;
            p2[i+1]=p2[i]*t%mod;
        }
        LL ans=0;
        rep(i,0,d+1) {
            LL t=g[i]*g[d-i]%mod*p1[i]%mod*p2[d-i]%mod*a[i]%mod;
            if ((d-i)&1) ans=(ans-t+mod)%mod;
            else ans=(ans+t)%mod;
        }
        return ans;
    }
    void init(int M) {///用到的最高次
        f[0]=f[1]=g[0]=g[1]=1;
        rep(i,2,M+5) f[i]=f[i-1]*i%mod;
        g[M+4]=powmod(f[M+4],mod-2);
        per(i,1,M+4) g[i]=g[i+1]*(i+1)%mod;///费马小定理筛逆元
    }

    ///函数用途：给出数列的（m+1）项，其中m为最高次方
    ///求出数列的前（n-1）项的和（从第0项开始）
    LL polysum(LL m,LL *a,LL n) { /// a[0].. a[m] \sum_{i=0}^{n-1} a[i]
        for(int i=0;i<=m;i++) b[i]=a[i];
        
        ///前n项和，其最高次幂加1
        b[m+1]=calcn(m,b,m+1);
        rep(i,1,m+2) b[i]=(b[i-1]+b[i])%mod;
        return calcn(m+1,b,n-1);
    }
    LL qpolysum(LL R,LL n,LL *a,LL m) { /// a[0].. a[m] \sum_{i=0}^{n-1} a[i]*R^i
        if (R==1) return polysum(n,a,m);
        a[m+1]=calcn(m,a,m+1);
        LL r=powmod(R,mod-2),p3=0,p4=0,c,ans;
        h[0][0]=0;h[0][1]=1;
        rep(i,1,m+2) {
            h[i][0]=(h[i-1][0]+a[i-1])*r%mod;
            h[i][1]=h[i-1][1]*r%mod;
        }
        rep(i,0,m+2) {
            LL t=g[i]*g[m+1-i]%mod;
            if (i&1) p3=((p3-h[i][0]*t)%mod+mod)%mod,p4=((p4-h[i][1]*t)%mod+mod)%mod;
            else p3=(p3+h[i][0]*t)%mod,p4=(p4+h[i][1]*t)%mod;
        }
        c=powmod(p4,mod-2)*(mod-p3)%mod;
        rep(i,0,m+2) h[i][0]=(h[i][0]+h[i][1]*c)%mod;
        rep(i,0,m+2) C[i]=h[i][0];
        ans=(calcn(m,C,n)*powmod(R,n)-c)%mod;
        if (ans<0) ans+=mod;
        return ans;
    }
}

 

本题的代码：

#include <bits/stdc++.h>

using namespace std;
typedef long long LL;

const LL mod=1e9+7;

/// 注意mod，使用前须调用一次 polysum::init(int M);
namespace polysum {
    #define rep(i,a,n) for (int i=a;i<n;i++)
    #define per(i,a,n) for (int i=n-1;i>=a;i--)
    const int D=1010000;///可能需要用到的最高次
    LL a[D],f[D],g[D],p[D],p1[D],p2[D],b[D],h[D][2],C[D];
    LL powmod(LL a,LL b){LL res=1;a%=mod;assert(b>=0);for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;}

     ///函数用途：给出数列的（d+1）项，其中d为最高次方项
    ///求出数列的第n项，数组下标从0开始
    LL calcn(int d,LL *a,LL n) { /// a[0].. a[d]  a[n]
        if (n<=d) return a[n];
        p1[0]=p2[0]=1;
        rep(i,0,d+1) {
            LL t=(n-i+mod)%mod;
            p1[i+1]=p1[i]*t%mod;
        }
        rep(i,0,d+1) {
            LL t=(n-d+i+mod)%mod;
            p2[i+1]=p2[i]*t%mod;
        }
        LL ans=0;
        rep(i,0,d+1) {
            LL t=g[i]*g[d-i]%mod*p1[i]%mod*p2[d-i]%mod*a[i]%mod;
            if ((d-i)&1) ans=(ans-t+mod)%mod;
            else ans=(ans+t)%mod;
        }
        return ans;
    }
    void init(int M) {///用到的最高次
        f[0]=f[1]=g[0]=g[1]=1;
        rep(i,2,M+5) f[i]=f[i-1]*i%mod;
        g[M+4]=powmod(f[M+4],mod-2);
        per(i,1,M+4) g[i]=g[i+1]*(i+1)%mod;///费马小定理筛逆元
    }

    ///函数用途：给出数列的（m+1）项，其中m为最高次方
    ///求出数列的前（n-1）项的和（从第0项开始）
    LL polysum(LL m,LL *a,LL n) { /// a[0].. a[m] \sum_{i=0}^{n-1} a[i]
        for(int i=0;i<=m;i++) b[i]=a[i];
        
        ///前n项和，其最高次幂加1
        b[m+1]=calcn(m,b,m+1);
        rep(i,1,m+2) b[i]=(b[i-1]+b[i])%mod;
        return calcn(m+1,b,n-1);
    }
    LL qpolysum(LL R,LL n,LL *a,LL m) { /// a[0].. a[m] \sum_{i=0}^{n-1} a[i]*R^i
        if (R==1) return polysum(n,a,m);
        a[m+1]=calcn(m,a,m+1);
        LL r=powmod(R,mod-2),p3=0,p4=0,c,ans;
        h[0][0]=0;h[0][1]=1;
        rep(i,1,m+2) {
            h[i][0]=(h[i-1][0]+a[i-1])*r%mod;
            h[i][1]=h[i-1][1]*r%mod;
        }
        rep(i,0,m+2) {
            LL t=g[i]*g[m+1-i]%mod;
            if (i&1) p3=((p3-h[i][0]*t)%mod+mod)%mod,p4=((p4-h[i][1]*t)%mod+mod)%mod;
            else p3=(p3+h[i][0]*t)%mod,p4=(p4+h[i][1]*t)%mod;
        }
        c=powmod(p4,mod-2)*(mod-p3)%mod;
        rep(i,0,m+2) h[i][0]=(h[i][0]+h[i][1]*c)%mod;
        rep(i,0,m+2) C[i]=h[i][0];
        ans=(calcn(m,C,n)*powmod(R,n)-c)%mod;
        if (ans<0) ans+=mod;
        return ans;
    }
}

LL num[1000010];
int main()
{

    LL n;
    int k;
    scanf("%lld%d",&n,&k);
    polysum::init(k+10);

    ///预处理前k+1个（从第0项开始），因为最高次为k+1次
    for(int i=0;i<=k+1;i++)
        num[i]=polysum::powmod((LL)i+1,k); ///注意此处是(i+1)^k

        ///计算前n项和
    LL ans=polysum::polysum(k+1,num,n)%mod;
    printf("%lld\n",ans);
    return 0;
}