[luoguP4705]玩游戏
前言
好的思路題
題目相關
傳送門
題意簡介
給定長度為nnn的序列aaa和長度為mmm的序列bbb
對所有k∈[1,t]k\in[1,t]k∈[1,t],求在兩個序列里分別隨機取一個數的和的kkk次方的期望,模998244353998244353998244353
數據范圍
1≤n,m≤105,1≤t≤1051\le n,m\le10^5,1\le t\le 10^51≤n,m≤105,1≤t≤105
題解
化一下式子,求期望相當于求所有情況的和取平均值
∑(a+b)k=∑∑i=0k(ki)aibk?i=∑i=0k(ki)(∑ai)(∑bk?i)=∑i=0k(ki)(∑ai)(∑bk?i)=∑i=0kk!i!(k?i)!(∑ai)(∑bk?i)=∑i=0kk!∑aii!∑bk?i(k?i)!=k!∑i=0k∑aii!∑bk?i(k?i)!\begin{aligned} \sum(a+b)^k&=\sum\sum_{i=0}^k\binom{k}{i}a^ib^{k-i}\\ &=\sum_{i=0}^k\binom{k}{i}(\sum a^i)(\sum b^{k-i})\\ &=\sum_{i=0}^k\binom{k}{i}(\sum a^i)(\sum b^{k-i})\\ &=\sum_{i=0}^k\frac{k!}{i!(k-i)!}(\sum a^i)(\sum b^{k-i})\\ &=\sum_{i=0}^kk!\frac{\sum a^i}{i!}\frac{\sum b^{k-i}}{(k-i)!}\\ &=k!\sum_{i=0}^k\frac{\sum a^i}{i!}\frac{\sum b^{k-i}}{(k-i)!}\\ \end{aligned}∑(a+b)k?=∑i=0∑k?(ik?)aibk?i=i=0∑k?(ik?)(∑ai)(∑bk?i)=i=0∑k?(ik?)(∑ai)(∑bk?i)=i=0∑k?i!(k?i)!k!?(∑ai)(∑bk?i)=i=0∑k?k!i!∑ai?(k?i)!∑bk?i?=k!i=0∑k?i!∑ai?(k?i)!∑bk?i??
我們如果求出所有∑ai\sum a^i∑ai和∑bi\sum b^i∑bi
考慮生成函數
Fi(x)=∑j=0∞aijxjF_i(x)=\sum_{j=0}^{\infty}a_i^jx^jFi?(x)=j=0∑∞?aij?xj
根據等比數列求和得到
Fi(x)=11?aixF_i(x)=\frac{1}{1-a_ix}Fi?(x)=1?ai?x1?
我們對所有生成函數求和
∑i=1nFi(x)=∑i=1n11?aix\sum_{i=1}^nF_i(x)=\sum_{i=1}^n\frac{1}{1-a_ix}i=1∑n?Fi?(x)=i=1∑n?1?ai?x1?
然后分治求分子、分母,用多項式求逆將分母乘進分子里
上面的式子算完后別忘了由于是期望所以要除以n*m
算法總復雜度O(nlog2n)\mathcal O(nlog^2n)O(nlog2n)
代碼
這份代碼優化了DFT與IDFT的次數,速度大大提高
#include<cstdio> #include<cctype> #include<cstring> #include<cstdlib> #include<vector> namespace fast_IO {const int IN_LEN=10000000,OUT_LEN=10000000;char ibuf[IN_LEN],obuf[OUT_LEN],*ih=ibuf+IN_LEN,*oh=obuf,*lastin=ibuf+IN_LEN,*lastout=obuf+OUT_LEN-1;inline char getchar_(){return (ih==lastin)&&(lastin=(ih=ibuf)+fread(ibuf,1,IN_LEN,stdin),ih==lastin)?EOF:*ih++;}inline void putchar_(const char x){if(oh==lastout)fwrite(obuf,1,oh-obuf,stdout),oh=obuf;*oh++=x;}inline void flush(){fwrite(obuf,1,oh-obuf,stdout);} } using namespace fast_IO; #define getchar() getchar_() #define putchar(x) putchar_((x)) typedef long long ll; #define rg register template <typename T> inline T max(const T a,const T b){return a>b?a:b;} template <typename T> inline T min(const T a,const T b){return a<b?a:b;} template <typename T> inline T mind(T&a,const T b){a=a<b?a:b;} template <typename T> inline T maxd(T&a,const T b){a=a>b?a:b;} template <typename T> inline T abs(const T a){return a>0?a:-a;} template <typename T> inline void swap(T&a,T&b){T c=a;a=b;b=c;} template <typename T> inline void swap(T*a,T*b){T c=a;a=b;b=c;} template <typename T> inline T gcd(const T a,const T b){if(!b)return a;return gcd(b,a%b);} template <typename T> inline T square(const T x){return x*x;}; template <typename T> inline void read(T&x) {char cu=getchar();x=0;bool fla=0;while(!isdigit(cu)){if(cu=='-')fla=1;cu=getchar();}while(isdigit(cu))x=x*10+cu-'0',cu=getchar();if(fla)x=-x; } template <typename T> void printe(const T x) {if(x>=10)printe(x/10);putchar(x%10+'0'); } template <typename T> inline void print(const T x) {if(x<0)putchar('-'),printe(-x);else printe(x); } const int maxn=524288,mod=998244353; inline int Md(const int x){return x>=mod?x-mod:x;} template<typename T> inline int pow(int x,T y) {rg int res=1;x%=mod;for(;y;y>>=1,x=(ll)x*x%mod)if(y&1)res=(ll)res*x%mod;return res; } namespace Poly///namespace of Poly { int W_[maxn],FW_[maxn],ha[maxn],hb[maxn],Inv[maxn]; inline void init(const int x) {rg int tim=0,lenth=1;while(lenth<x)lenth<<=1,tim++;for(rg int i=1;i<lenth;i++){W_[i]=pow(3,(mod-1)/i/2);FW_[i]=pow(W_[i],mod-2);} // Inv[0]=Inv[1]=1;for(rg int i=2;i<x;i++)Inv[i]=(ll)(mod-mod/i)*Inv[mod%i]%mod; } int L; inline void NTT(int*A,const int fla)//prepare:init L {for(rg int i=0,j=0;i<L;i++){if(i>j)swap(A[i],A[j]);for(rg int k=L>>1;(j^=k)<k;k>>=1);}for(rg int i=1;i<L;i<<=1){const int w=fla==-1?FW_[i]:W_[i];for(rg int j=0,J=i<<1;j<L;j+=J){int K=1;for(rg int k=0;k<i;k++,K=(ll)K*w%mod){const int x=A[j+k],y=(ll)A[j+k+i]*K%mod;A[j+k]=Md(x+y),A[j+k+i]=Md(mod+x-y);}}} } inline int Quadratic_residue(const int a) {if(a==0)return 0;int b=(rand()<<14^rand())%mod;while(pow(b,(mod-1)>>1)!=mod-1)b=(rand()<<14^rand())%mod;int s=mod-1,t=0,x,inv=pow(a,mod-2),f=1;while(!(s&1))s>>=1,t++,f<<=1;t--,x=pow(a,(s+1)>>1),f>>=1;while(t){f>>=1;if(pow((int)((ll)inv*x%mod*x%mod),f)!=1)x=(ll)x*pow(b,s)%mod;t--,s<<=1;}return min(x,mod-x); } struct poly {std::vector<int>A;poly(){A.resize(0);}poly(const int x){A.resize(1),A[0]=x;}inline int&operator[](const int x){return A[x];}inline int operator[](const int x)const{return A[x];}inline void clear(){A.clear();}inline unsigned int size()const{return A.size();}inline void resize(const unsigned int x){A.resize(x);}void RE(const int x){A.resize(x);for(rg int i=0;i<x;i++)A[i]=0; }void readin(const int MAX){A.resize(MAX);for(rg int i=0;i<MAX;i++)read(A[i]);}void putout()const{for(rg unsigned int i=0;i<A.size();i++)print(A[i]),putchar(' ');}inline poly operator +(const poly b)const{poly RES;RES.resize(max(size(),b.size()));for(rg unsigned int i=0;i<RES.size();i++)RES[i]=Md((i<size()?A[i]:0)+(i<b.size()?b[i]:0));return RES;}inline poly operator -(const poly b)const{poly RES;RES.resize(max(size(),b.size()));for(rg unsigned int i=0;i<RES.size();i++)RES[i]=Md((i<size()?A[i]:0)+mod-(i<b.size()?b[i]:0));return RES;}inline poly operator *(const int b)const{poly RES=*this;for(rg unsigned int i=0;i<RES.size();i++)RES[i]=(ll)RES[i]*b%mod;return RES;}inline poly operator /(const int b)const{poly RES=(*this)*pow(b,mod-2);return RES;}inline poly operator *(const poly b)const{const int RES=A.size()+b.size()+1;L=1;while(L<RES)L<<=1;poly c;c.A.resize(RES);memset(ha,0,sizeof(int)*L);memset(hb,0,sizeof(int)*L);for(rg unsigned int i=0;i<A.size();i++)ha[i]=A[i];for(rg unsigned int i=0;i<b.A.size();i++)hb[i]=b.A[i];NTT(ha,1),NTT(hb,1);for(rg int i=0;i<L;i++)ha[i]=(ll)ha[i]*hb[i]%mod;NTT(ha,-1);const int inv=pow(L,mod-2);for(rg int i=0;i<RES;i++)c.A[i]=(ll)ha[i]*inv%mod;return c;}inline poly inv()const{poly C;if(A.size()==1){C=*this;C[0]=pow(C[0],mod-2);return C;}C.resize((A.size()+1)>>1);for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i];C=C.inv();L=1;while(L<(int)size()*2-1)L<<=1;for(rg unsigned int i=0;i<A.size();i++)ha[i]=A[i];for(rg unsigned int i=0;i<C.size();i++)hb[i]=C[i];memset(ha+A.size(),0,sizeof(int)*(L-A.size()));memset(hb+C.size(),0,sizeof(int)*(L-C.size()));NTT(ha,1),NTT(hb,1);for(rg int i=0;i<L;i++)ha[i]=(ll)(2-(ll)hb[i]*ha[i]%mod+mod)*hb[i]%mod;NTT(ha,-1);const int inv=pow(L,mod-2);C.resize(size());for(rg unsigned int i=0;i<size();i++)C[i]=(ll)ha[i]*inv%mod;return C;} /* inline poly inv()const{poly C;if(A.size()==1){C=*this;C[0]=pow(C[0],mod-2);return C;}C.resize((A.size()+1)>>1);for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i];C=C.inv();poly D=(poly)2-*this*C;D.resize(size());D=D*C;D.resize(size());return D;}*///大常數版本 inline void Reverse(const int n){A.resize(n);for(rg int i=0,j=n-1;i<j;i++,j--)swap(A[i],A[j]);}inline poly operator /(const poly B)const{poly a=*this,b=B;a.Reverse(size()),b.Reverse(B.size());b.resize(size()-B.size()+1);b=b.inv();b=b*a;b.Reverse(size()-B.size()+1);return b;}inline poly operator %(const poly B)const{poly C=(*this)-(*this)/B*B;C.resize(B.size()-1);return C;}inline poly diff()const{poly C;C.resize(size()-1);for(rg unsigned int i=1;i<size();i++)C[i-1]=(ll)A[i]*i%mod;return C;}inline poly inte()const{poly C;C.resize(size()+1);for(rg unsigned int i=0;i<size();i++)C[i+1]=(ll)A[i]*Inv[i+1]%mod;C[0]=0;return C;}inline poly ln()const{poly C=(diff()*inv()).inte();C.resize(size());return C;}inline poly exp()const{poly C;if(size()==1){C=*this;C[0]=1;return C;}C.resize((size()+1)>>1);for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i];C=C.exp();C.resize(size());poly D=(poly)1-C.ln()+*this;D=D*C;D.resize(size());return D;}inline poly sqrt()const{poly C;if(size()==1){C=*this;C[0]=Quadratic_residue(C[0]);return C;}C.resize((size()+1)>>1);for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i];C=C.sqrt();C.resize(size());C=(C+*this*C.inv())*(int)499122177;C.resize(size());return C;}inline poly operator >>(const unsigned int x)const{poly C;if(size()<x){C.resize(0);return C;}C.resize(size()-x);for(rg unsigned int i=0;i<C.size();i++)C[i]=A[i+x];return C;}inline poly operator <<(const unsigned int x)const{poly C;C.RE(size()+x);for(rg unsigned int i=0;i<size();i++)C[i+x]=A[i];return C;}inline poly Pow(const unsigned int x)const{for(rg unsigned int i=0;i<size();i++)if(A[i]){poly C=((((*this/A[i])>>i).ln()*x).exp()*pow(A[i],x))<<(min(i*x,size()));C.resize(size());return C;}return *this;}inline void cheng(const poly&B){for(rg unsigned int i=0;i<size();i++)A[i]=(ll)A[i]*B[i]%mod; }inline void jia(const poly&B){for(rg unsigned int i=0;i<size();i++)A[i]=Md(A[i]+B[i]); }inline void DFT(){memset(ha,0,sizeof(int)*L);for(rg unsigned int i=0;i<A.size();i++)ha[i]=A[i];NTT(ha,1);resize(L);for(rg int i=0;i<L;i++)A[i]=ha[i];}inline void IDFT(){memset(ha,0,sizeof(int)*L);for(rg unsigned int i=0;i<A.size();i++)ha[i]=A[i];NTT(ha,-1);const int inv=pow(L,mod-2);for(rg int i=0;i<L;i++)A[i]=(ll)ha[i]*inv%mod;while(size()&&!A[size()-1])A.pop_back();} }; }///namespace of Poly Poly::poly a[maxn],za[maxn],b[maxn],zb[maxn],A,B; int n,m,t; void fz(Poly::poly*mu,Poly::poly*zi,const int l,const int r) {if(l==r)return;const int mid=(l+r)>>1;fz(mu,zi,l,mid),fz(mu,zi,mid+1,r);const int RES=mu[l].size()<<1|1;Poly::L=1;while(Poly::L<RES)Poly::L<<=1;const int ll=l,rr=mid+1;zi[ll].DFT(),mu[ll].DFT(),zi[rr].DFT(),mu[rr].DFT();zi[ll].cheng(mu[rr]);zi[rr].cheng(mu[ll]);zi[ll].jia(zi[rr]);zi[ll].IDFT();mu[ll].cheng(mu[rr]);mu[ll].IDFT(); } int fac[maxn+1],inv[maxn+1]; int main() {fac[0]=1;for(rg int i=1;i<=maxn;i++)fac[i]=(ll)fac[i-1]*i%mod;inv[maxn]=pow(fac[maxn],mod-2);for(rg int i=maxn;i>=1;i--)inv[i-1]=(ll)inv[i]*i%mod;Poly::init(maxn);///namespace of Polyread(n),read(m);for(rg int i=1;i<=n;i++)a[i].resize(2),a[i][0]=1,read(a[i][1]),a[i][1]=Md(mod-a[i][1]),za[i].resize(1),za[i][0]=1;for(rg int i=1;i<=m;i++)b[i].resize(2),b[i][0]=1,read(b[i][1]),b[i][1]=Md(mod-b[i][1]),zb[i].resize(1),zb[i][0]=1;for(rg int i=1;i<=m;i++);read(t);fz(a,za,1,n);fz(b,zb,1,m);a[1].resize(100001);b[1].resize(100001);A=za[1]*a[1].inv();B=zb[1]*b[1].inv();A.resize(100001);B.resize(100001);for(rg int i=1;i<=100000;i++)A[i]=(ll)A[i]*inv[i]%mod,B[i]=(ll)B[i]*inv[i]%mod;A=A*B;const int INV=pow((ll)n*m%mod,mod-2);for(rg int i=1;i<=t;i++)print((ll)A[i]*fac[i]%mod*INV%mod),putchar('\n');return flush(),0; }題解2
這題有更清真的做法
求所有∑ai\sum a^i∑ai和∑bi\sum b^i∑bi的時候能用ln優化
考慮生成函數
Fi(x)=∑j=0∞aijxjF_i(x)=\sum_{j=0}^{\infty}a_i^jx^jFi?(x)=j=0∑∞?aij?xj
根據等比數列求和得到
Fi(x)=11?aixF_i(x)=\frac{1}{1-a_ix}Fi?(x)=1?ai?x1?
我們對所有生成函數求和,設
f(x)=∑i=1nFi(x)=∑i=1n11?aixf(x)=\sum_{i=1}^nF_i(x)=\sum_{i=1}^n\frac{1}{1-a_ix}f(x)=i=1∑n?Fi?(x)=i=1∑n?1?ai?x1?
直接求常數很大
我們考慮求導
眾所周知(ln(x))′=1x,(ax)′=a(ln(x))'=\frac1x,(ax)'=a(ln(x))′=x1?,(ax)′=a
根據復合函數求導法則
(ln(1?aix))′=?ai1?aix(ln(1-a_ix))'=\frac{-a_i}{1-a_ix}(ln(1?ai?x))′=1?ai?x?ai??
設g(x)=∑i=1n?ai1?aixg(x)=\sum_{i=1}^n\frac{-a_i}{1-a_ix}g(x)=i=1∑n?1?ai?x?ai??
我們發現f(x)+g(x)x=nf(x)+g(x)x=nf(x)+g(x)x=n
移項得f(x)=n?g(x)xf(x)=n-g(x)xf(x)=n?g(x)x
那么我們只要求出g(x)g(x)g(x)就能求出f(x)f(x)f(x)
g(x)=∑i=1n?ai1?aix=∑i=1n(ln(1?aix))′=(∑i=1nln(1?aix))′=(ln(∏i=1n(1?aix)))′\begin{aligned} g(x)&=\sum_{i=1}^n\frac{-a_i}{1-a_ix}\\ &=\sum_{i=1}^n(ln(1-a_ix))'\\ &=(\sum_{i=1}^nln(1-a_ix))'\\ &=(ln(\prod_{i=1}^n(1-a_ix)))'\\ \end{aligned}g(x)?=i=1∑n?1?ai?x?ai??=i=1∑n?(ln(1?ai?x))′=(i=1∑n?ln(1?ai?x))′=(ln(i=1∏n?(1?ai?x)))′?
這個算法總復雜度也是O(nlog2n)\mathcal O(nlog^2n)O(nlog2n)
不貼代碼了,我寫的ln常數太大了,過不了
懶得卡常了
總結
一道經典模型好題
總結
以上是生活随笔為你收集整理的[luoguP4705]玩游戏的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 指数型生成函数[bzoj3456]城市规
- 下一篇: 多项式快速插值