[bzoj3157][bzoj3516][bzoj4126]国王奇遇记

Description

求

$$\sum_{i=1}^{n}i^mm^i$$
n<=10^9,普通版m<=200，加强版m<=1000，加强版之再加强版m<=500000

Solution

1.普通版

很容易想到倍增的思路，只不过很难实现。
设$F(n,k)=\sum_{i=1}^{n}i^km^i$
然后,

$$F(n+1,k)=\sum_{i=1}^{n+1}i^km^i$$
$$=m+\sum_{i=2}^{n+1}i^km^i$$
$$=m+\sum_{i=1}^{n}(i+1)^km^{i+1}$$
$$=m+m\sum_{i=1}^{n}m^i\sum_{j=0}^{k}C_k^ji^j$$
$$=m+m\sum_{j=0}^{k}C_k^j\sum_{i=1}^{n}i^jm^i$$
$$=m+m\sum_{j=0}^{k}C_k^jF(n,j)$$
且 $$F(2n,k)=\sum_{i=1}^{2n}i^km^i$$
$$=\sum_{i=1}^{n}i^km^i+\sum_{i=n+1}^{2n}i^km^i$$
$$=F(n,k)+\sum_{i=1}^{n}(i+n)^km^{i+n}$$
$$=F(n,k)+m^n\sum_{i=1}^{n}m^i\sum_{j=0}^{k}C_k^ji^jn^{k-j}$$
$$=F(n,k)+m^n\sum_{j=0}^{k}C_k^jn^{k-j}\sum_{i=1}^{n}i^jm^i$$
$$=F(n,k)+m^n\sum_{j=0}^{k}C_k^jn^{k-j}F(n,j)$$
就这样每次m^2转移就好了。

Code

#include<map>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define N 205
using namespace std;
typedef long long ll;
const int mo=1000000000+7;
ll f[52][N],c[N][N],n;
int m;
ll mi(ll x,ll y) {
    ll z=x;
    for(y--;y;y/=2) {
        if (y&1) (z*=x)%=mo;
        (x*=x)%=mo;
    }
    return z;
}
void calc(ll x,int y) {
    if (x==1) {
        fo(i,0,m) f[y][i]=m;
        return;
    }
    if (x&1) {
        calc(x-1,y+1);
        fo(i,0,m) {
            fo(j,0,i) (f[y][i]+=c[i][j]*f[y+1][j])%=mo;
            (f[y][i]*=m)%=mo;(f[y][i]+=m)%=mo;
        }
    } else {
        calc(x/=2,y+1);ll pow=mi(m,x);
        fo(i,0,m) {
            ll sum=1;
            fo(j,0,i) (f[y][i]+=sum*c[i][i-j]%mo*f[y+1][i-j])%=mo,(sum*=x)%=mo;
            (f[y][i]*=pow)%=mo;(f[y][i]+=f[y+1][i])%=mo;
        }
    }
}
int main() {
    scanf("%lld%d",&n,&m);c[0][0]=1;
    fo(i,1,m) {
        c[i][0]=1;
        fo(j,1,i) (c[i][j]=c[i-1][j]+c[i-1][j-1])%=mo;
    } 
    calc(n,1);
    printf("%lld",f[1][m]);
}

2.加强版

奥妙重重的开头，直接设$S(k)=\sum_{i=1}^{n}i^km^i$
那么

$$mS(k)-S(k)=\sum_{i=1}^{n}i^km^{i+1}-\sum_{i=1}^{n}i^km^i$$
$$=\sum_{i=2}^{n+1}(i-1)^km^i-\sum_{i=1}^{n}i^km^i$$
$$=n^km^{n+1}+\sum_{i=1}^{n}m^i((i-1)^k-i^k)$$
$$=n^km^{n+1}+\sum_{i=1}^{n}m^i\sum_{j=0}^{k-1}C_k^ji^j(-1)^{k-j}$$
$$=n^km^{n+1}+\sum_{j=0}^{k-1}(-1)^{k-j}C_k^j\sum_{i=1}^{n}i^jm^i$$
$$=n^km^{n+1}+\sum_{j=0}^{k-1}(-1)^{k-j}C_k^jS(j)$$
直接递推即可。

Code

#include<cstdio>
#include<cstring>
#include<algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define N 1005
using namespace std;
typedef long long ll;
const int mo=1000000000+7;
ll f[N],c[N][N],n,m,sum,ni;
ll mi(ll x,ll y) {
    ll z=x;
    for(y--;y;y/=2) {
        if (y&1) (z*=x)%=mo;
        (x*=x)%=mo;
    }
    return z;
}
int main() {
    scanf("%lld%lld",&n,&m);
    if (m==1) {
        printf("%lld",n*(n+1)/2%mo);
        return 0;
    }
    sum=mi(m,n+1);ni=mi(m-1,mo-2);
    c[0][0]=1;
    fo(i,1,m) {
        c[i][0]=1;
        fo(j,1,i) c[i][j]=(c[i-1][j]+c[i-1][j-1])%mo;
    }
    f[0]=(sum-m+mo)%mo;(f[0]*=ni)%=mo;
    fo(i,1,m) {
        (sum*=n)%=mo;f[i]=sum;
        fo(j,0,i-1) {
            ll z=((i-j)&1)?-1:1;
            (f[i]+=z*c[i][j]*f[j]%mo)%=mo;
        }
        (f[i]+=mo)%=mo;(f[i]*=ni)%=mo;
    }
    printf("%lld",f[m]);
}

3.加强版之再加强版

表示一脸懵逼，不会做，在这里贴一个别人的证明题解
%%%（下面的超链接请戳这里）
好像很厉害的样子呀！！！！

Code

#include<cstdio>
typedef long long ll;
#define N 500010
#define p 1000000007
int n,m,im,fac[N],i,j,inv[N],a[N],b[N],fa,fb,f0,pr[N],pt,pwm[N],vis[N],flag,tmp,ans;
int power(int t,int k){
    int f=1;
    for(;k;k>>=1,t=1LL*t*t%p)if(k&1)f=1LL*f*t%p;
    return f;
}
#define C(n,k) (1LL*fac[n]*inv[k]%p*inv[n-k]%p)
int main(){
    scanf("%d%d",&n,&m);n++;
    if(m==1)return printf("%lld\n",(1LL*n*(n-1)/2)%p),0;
    for(pwm[1]=inv[1]=1,i=2;i<=m+1;i++){
        if(!vis[i])pr[++pt]=i,pwm[i]=power(i,m),inv[i]=power(i,p-2);
        for(j=1;j<=pt&&1LL*i*pr[j]<=m+1;j++){
            vis[i*pr[j]]=1,
            pwm[i*pr[j]]=1LL*pwm[i]*pwm[pr[j]]%p,
            inv[i*pr[j]]=1LL*inv[i]*inv[pr[j]]%p;
            if(i%pr[j]==0)break;
        }
    }
    for(fac[0]=inv[0]=i=1;i<=m+1;i++)fac[i]=1LL*fac[i-1]*i%p,inv[i]=1LL*inv[i-1]*inv[i]%p;
    for(im=power(m,p-2),a[i=0]=1;i<=m;i++)a[i+1]=1LL*a[i]*im%p,b[i+1]=1LL*(pwm[i]+b[i])*im%p;
    for(flag=1,i=0;i<=m+1;i++,flag=-flag){
        tmp=C(m+1,i)*flag;if(tmp<0)tmp+=p;
        fa=(fa+1LL*tmp*a[i])%p,
        fb=(fb+1LL*tmp*b[i])%p;
    }
    f0=(p-1LL*fb*power(fa,p-2)%p)%p;
    if(n<=m){
        ans=(1LL*power(m,n)*(1LL*a[n]*f0%p+b[n])%p-f0+p)%p;
        printf("%d\n",ans);
    }
    else{
        for(flag=1,i=m;i>=0;i--,flag=-flag){
            tmp=(1LL*a[i]*f0+b[i])%p*inv[i]%p*inv[m-i]%p*power(n-i,p-2)%p;
            if(flag<0)tmp=p-tmp;ans=(ans+tmp)%p;
        }
        for(i=n-m;i<=n;i++)ans=1LL*ans*i%p;
        ans=(1LL*power(m,n)*ans%p-f0+p)%p;
        printf("%d\n",ans);
    }
}