【BZOJ4126】【BZOJ3516】【BZOJ3157】国王奇遇记 线性插值

题目描述

　　三倍经验题。

　　给你\(n,m\)，求
\[ \sum_{i=1}^ni^mm^i \]
　　\(n\leq {10}^9,1\leq m\leq 500000\)

题解

　　当\(m=1\)时\(ans=\frac{n(n+1)}{2}\)

　　剩下的部分这篇博客有讲YWW's Blog

　　时间复杂度：\(O(m+\log n)\)

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const ll p=1000000007;
ll fp(ll a,ll b)
{
    ll s=1;
    for(;b;b>>=1,a=a*a%p)
        if(b&1)
            s=s*a%p;
    return s;
}
int pri[100010];
int b[1000010];
int cnt;
ll s[1000010];
ll fac[1000010];
ll ifac[1000010];
ll inv[1000010];
ll f1[1000010];
ll f2[1000010];
ll f[1000010];
ll pre[1000010];
ll suf[1000010];
ll getc(int x,int y)
{
    return fac[x]*ifac[y]%p*ifac[x-y]%p;
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("bzoj4126.in","r",stdin);
    freopen("bzoj4126.out","w",stdout);
#endif
    int n,m;
    scanf("%d%d",&n,&m);
    if(m==1)
    {
        printf("%lld\n",ll(n)*(n+1)/2%p);
        return 0;
    }
    n++;
    fac[0]=fac[1]=ifac[0]=ifac[1]=inv[1]=1;
    for(int i=2;i<=m+2;i++)
    {
        inv[i]=-p/i*inv[p%i]%p;
        ifac[i]=ifac[i-1]*inv[i]%p;
        fac[i]=fac[i-1]*i%p;
    }
    s[0]=0;
    s[1]=1;
    for(int i=2;i<=m+2;i++)
    {
        if(!b[i])
        {
            s[i]=fp(i,m);
            pri[++cnt]=i;
        }
        for(int j=1;j<=cnt&&i*pri[j]<=m+2;j++)
        {
            b[i*pri[j]]=1;
            s[i*pri[j]]=s[i]*s[pri[j]]%p;
            if(i%pri[j]==0)
                break;
        }
    }
    f1[0]=1;
    f2[0]=0;
    ll invm=fp(m,p-2);
    for(int i=1;i<=m+1;i++)
    {
        f1[i]=f1[i-1]*invm%p;
        f2[i]=(f2[i-1]+s[i-1])*invm%p;
    }
    ll v1=0,v2=0;
    for(int i=0;i<=m+1;i++)
    {
        v1=(v1+((m+1-i)&1?-1:1)*getc(m+1,i)*f1[i])%p;
        v2=(v2+((m+1-i)&1?-1:1)*getc(m+1,i)*f2[i])%p;
    }
    f[0]=-v2*fp(v1,p-2)%p;
    for(int i=1;i<=m+1;i++)
        f[i]=(f1[i]*f[0]+f2[i])%p;
    if(n<=m+1)
    {
        ll ans=fp(m,n)*f[n]-f[0];
        ans=(ans%p+p)%p;
        printf("%lld\n",ans);
        return 0;
    }
    for(int i=0;i<=m;i++)
    {
        pre[i]=n-i;
        if(i)
            pre[i]=pre[i-1]*pre[i]%p;
    }
    for(int i=m;i>=0;i--)
    {
        suf[i]=n-i;
        if(i!=m)
            suf[i]=suf[i+1]*suf[i]%p;
    }
    ll ans=0;
    for(int i=0;i<=m;i++)
    {
        ll v=1;
        if(i)
            v=v*pre[i-1]%p;
        if(i!=m)
            v=v*suf[i+1]%p;
        ans=(ans+f[i]*v%p*ifac[i]%p*ifac[m-i]%p*((m-i)&1?-1:1))%p;
    }
    ans=fp(m,n)*ans-f[0];
    ans=(ans%p+p)%p;
    printf("%lld\n",ans);
    return 0;
}

转载于:https://www.cnblogs.com/ywwyww/p/8600445.html