简要题意:
你有一个数,初始为 0 0 0,有 n n n 个机会,每个机会有 Q / ( 1 + Q ) Q/(1+Q) Q/(1+Q) 的概率使你的数 +1,请计算所有小于等于你的数的自然数的 k k k 次幂之和的期望。
题解:
设 f n , i f_{n,i} fn,i 表示用了前 n n n 个机会,你的数的 i i i 次幂的期望。
转移考虑二项式展开,乘的东西不变,直接多项式快速幂。
根据期望的线性性,由于 k k k 次幂之和是个 k + 1 k+1 k+1 次多项式,插值求出它的系数表达然后线性相加即可。
代码:
#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const
using std::cerr;
using std::cout;
#define fi first
#define se second
cs int mod=1e9+7;
inline int add(int a,int b){return a+b>=mod?a+b-mod:a+b;}
inline int dec(int a,int b){return a-b<0?a-b+mod:a-b;}
inline int mul(int a,int b){ll r=(ll)a*b;return r>=mod?r%mod:r;}
inline void Inc(int &a,int b){a+=b-mod;a+=a>>31&mod;}
inline void Dec(int &a,int b){a-=b;a+=a>>31&mod;}
inline void Mul(int &a,int b){a=mul(a,b);}
inline int po(int a,int b){int r=1;for(;b;b>>=1,Mul(a,a))if(b&1)Mul(r,a);return r;}
inline void ex_gcd(int a,int b,int&x,int &y){
if(!b){x=1,y=0;return;}ex_gcd(b,a%b,y,x);y-=a/b*x;
}inline int Inv(int a){int x,y;ex_gcd(mod,a,y,x);return x+(x>>31&mod);}
cs int bit=19,SIZE=1<<bit|7;
int fac[SIZE],ifc[SIZE],inv[SIZE];
void init_fac(){
fac[0]=fac[1]=1;
ifc[0]=ifc[1]=1;
inv[0]=inv[1]=1;
for(int re i=2;i<SIZE;++i){
fac[i]=mul(fac[i-1],i);
inv[i]=mul(inv[mod%i],mod-mod/i);
ifc[i]=mul(ifc[i-1],inv[i]);
}
}
namespace MTT{
cs double PI=acos(-1),PI2=PI+PI;
struct cp{
double x,y;cp(){}
cp(double _x,double _y):x(_x),y(_y){}
friend cp operator+(cs cp &a,cs cp &b){return cp(a.x+b.x,a.y+b.y);}
friend cp operator-(cs cp &a,cs cp &b){return cp(a.x-b.x,a.y-b.y);}
friend cp operator*(cs cp &a,cs cp &b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
cp& operator+=(cs cp &b){x+=b.x,y+=b.y;return *this;}
cp& operator-=(cs cp &b){x-=b.x,y-=b.y;return *this;}
cp& operator*=(cs cp &b){return *this=*this*b;}
cp conj()cs{return cp(x,-y);}
};
inline cp omega(int i,int k){
return cp(cos(PI2*i/k),sin(PI2*i/k));
}
int r[SIZE];cp *w[bit+1];
void init_omega(){
for(int re i=1;i<=bit;++i)
w[i]=new cp[1<<(i-1)];
for(int re d=1;d<=bit;++d){
cp wn=omega(1,1<<d);
for(int re i=0;i<(1<<(d-1));++i)
w[d][i]=(i&31)?w[d][i-1]*wn:omega(i,1<<d);
}
}
void DFT(cp *A,int l){
for(int re i=0;i<l;++i)
if(i<r[i])std::swap(A[i],A[r[i]]);
for(int re i=1,d=1;i<l;i<<=1,++d)
for(int re j=0;j<l;j+=i<<1)
for(int re k=0;k<i;++k){
cp &t1=A[j+k],&t2=A[i+j+k];
cp t=t2*w[d][k];t2=t1-t,t1+=t;
}
}
void init_len(int l){
for(int re i=1;i<l;++i)
r[i]=r[i>>1]>>1|((i&1)?l>>1:0);
}
void mul(int *a,int *b,int l,int *c){
static cp A[SIZE],B[SIZE],C[SIZE],D[SIZE];
for(int re i=0;i<l;++i){
A[i]=cp(a[i]&0x7fff,a[i]>>15);
B[i]=cp(b[i]&0x7fff,b[i]>>15);
}init_len(l),DFT(A,l),DFT(B,l);
for(int re i=0;i<l;++i){
int u=(l-i)&(l-1);
C[i]=(A[i].conj()+A[u])*cp(.5,0)*B[u];
D[i]=(A[i].conj()-A[u])*cp(0,.5)*B[u];
}DFT(C,l);DFT(D,l);
for(int re i=0;i<l;++i){
ll x=C[i].x/l+.5,y=(C[i].y+D[i].x)/l+.5,z=D[i].y/l+.5;
x%=mod,y%=mod,z%=mod;c[i]=(x+(y<<15)%mod+(z<<30)%mod)%mod;
}
}
}
typedef std::vector<int> Poly;
Poly operator*(Poly a,int b){
for(size_t re i=0;i<a.size();++i)
Mul(a[i],b);return a;
}Poly operator+(cs Poly &a,cs Poly &b){
Poly c=a;if(c.size()<b.size())c.resize(b.size());
for(size_t re i=0;i<b.size();++i)Inc(c[i],b[i]);
return c;
}Poly operator-(cs Poly &a,cs Poly &b){
Poly c=a;if(c.size()<b.size())c.resize(b.size());
for(size_t re i=0;i<b.size();++i)Dec(c[i],b[i]);
return c;
}Poly operator*(Poly a,Poly b){
if(!a.size()||!b.size())
return Poly(0);
int deg=a.size()+b.size()-1;
int l=1;while(l<deg)l<<=1;
a.resize(l),b.resize(l);
MTT::mul(&a[0],&b[0],l,&a[0]);
return Poly(a.begin(),a.begin()+deg);
}Poly Inv(cs Poly &a,int lim){
int n=a.size();Poly c,b(1,Inv(a[0]));
for(int re l=4;(l>>2)<lim;l<<=1){
c.resize(l>>1);
for(int re i=0;i<(l>>1);++i)
c[i]=i<n?a[i]:0;
b=b*2-c*b*b;b.resize(l>>1);
}return Poly(b.begin(),b.begin()+lim);
}Poly Integ(Poly a){
if(!a.size())return Poly(0);
a.push_back(0);
for(size_t re i=a.size()-1;i;--i)
a[i]=mul(a[i-1],inv[i]);
a[0]=0;return a;
}Poly Deriv(Poly a){
if(!a.size())return Poly(0);
for(size_t re i=1;i<a.size();++i)
a[i-1]=mul(a[i],i);
a.pop_back();return a;
}Poly Ln(Poly a,int lim){
a=Deriv(a)*Inv(a,lim);
a.resize(lim+1);return Integ(a);
}Poly Exp(cs Poly &a,int lim){
int n=a.size();Poly c,b(1,1);
for(int re i=2;(i>>1)<lim;i<<=1){
c=Ln(b,i);Dec(c[0],1);
for(int re j=0;j<i;++j)
c[j]=dec(j<n?a[j]:0,c[j]);
b=b*c;b.resize(i);
}return Poly(b.begin(),b.begin()+lim);
}
std::pair<Poly,Poly> interpolation(int l,int r,cs Poly &A,cs Poly &coef){
if(l==r){
Poly h(2,0);h[0]=dec(0,l),h[1]=1;
Poly g(1,mul(A[l],coef[l]));
return std::pair<Poly,Poly>(h,g);
}int mid=(l+r)>>1;
std::pair<Poly,Poly> a=interpolation(l,mid,A,coef);
std::pair<Poly,Poly> b=interpolation(mid+1,r,A,coef);
return std::pair<Poly,Poly>(a.fi*b.fi,a.fi*b.se+b.fi*a.se);
}Poly interpolation(cs Poly &G){
int n=G.size();Poly coef(n);
for(int re i=0;i<n;++i)
coef[i]=mul(ifc[i],((n-i)&1)?ifc[n-i-1]:mod-ifc[n-i-1]);
return interpolation(0,n-1,G,coef).se;
}
Poly calc(int k){//f(x)=\sum_{i=0}^x i^k
Poly f(k+2);
for(int i=0,sm=0;i<=k+1;++i)
f[i]=sm,Inc(sm,po(i+1,k));
return interpolation(f);
}
int n,k,Q;
void Main(){
MTT::init_omega();init_fac();
scanf("%d%d%d",&n,&k,&Q);Poly f(k+2);int iv=Inv(Q+1);
for(int re i=0;i<=k+1;++i)f[i]=mul(mul(Q,iv),ifc[i]);
Inc(f[0],iv);f=Ln(f,k+2);int vl=po(Q+1,n);
for(int re i=0;i<=k+1;++i)Mul(f[i],n);
f=Exp(f,k+2);Poly g=calc(k);int ans=0;
for(int re i=0;i<=k+1;++i)
Inc(ans,mul(mul(mul(f[i],vl),fac[i]),g[i]));
cout<<ans<<"\n";
}
inline void file(){
#ifdef zxyoi
freopen("qin.in","r",stdin);
#else
freopen("qin.in","r",stdin);
freopen("qin.out","w",stdout);
#endif
}signed main(){file();Main();return 0;}