前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >YbtOJ 976「母函数」随机减法

YbtOJ 976「母函数」随机减法

作者头像
yzxoi
发布2022-09-19 14:03:20
2.1K0
发布2022-09-19 14:03:20
举报
文章被收录于专栏:OI

YbtOJ 976「母函数」随机减法

题目链接:YbtOJ #976

小 A 有一个长度为 n 的序列 a 和一个初始值为 0 的计数器 cnt,他想要对其进行 k 次操作。

每次操作,他会等概率随机选中一个 i,将 a_i1,并将 cnt 加上 此时 除 a_i 以外所有数的乘积,即 \prod_{j\not=i}a_j

现在,他希望知道 cnt 在模 10^9+7 意义下的期望值。

1\le n\le5\times10^31\le k\le10^90\le a_i < 10^9+7

Solution

容易发现将 a_i1 后,除它以外所有数的乘积恰好是 \prod_{i=1}^na_i 的变化量。

所以说,答案实际上就是原本的 \prod_{i=1}^na_i 减去修改后 \prod_{i=1}^na_i’ 的期望值。

f_{i,j} 表示前 i 个数一共修改了 j 次的所有方案下乘积之和,则:

f_{i,p+q}=\sum C_{p+q}^p\times (a_i-p)\times f_{i-1,q}

于是:

\frac{f_{i,p+q}}{(p+q)!}=\sum\frac{a_i-p}{p!}\times\frac{f_{i-1,q}}{q!}

F_i(x)=\sum_{p=0}^{+\infty}f_{i,p}\frac{x^p}{p!},G_i(x)=\sum_{p=0}^{+\infty}(a_i-p)\frac{x^p}{p!},得到:

F_i(x)=F_{i-1}(x)*G_i(x)

因此只要把 G_{1\sim n}(x)n 个生成函数卷起来就能得到 F_n(x),而它的 k 次项系数就是 \frac{f_{n,k}}{k!} 了。

对于 G_i(x),我们把 a_i-p 分开来:

G_i(x)=a_i\sum_{p=0}^{+\infty}\frac{x^p}{p!}-\sum_{p=0}^{+\infty}\frac{x^{p+1}}{p!}=(a_i-x)e^x

A_i(x)=a_i-x,发现 F_n(x) 就是 A_{1\sim n}(x)n 个生成函数卷起来之后再卷上 e^{nx}

很容易 O(n^2) 暴力求出 A_{1\sim n}(x) 卷起来后每一项的系数 f_i,于是:

乘上一个 k! 得到了总和 f_{n,k},再除以总方案数 n^k 得到期望:

E=\sum_{i=0}^{k}\frac{f_i\times k^{\underline i}}{n^i}

最终答案就是\prod_{i=1}^na_i-E

注意一开始卷积直接暴力卷即可。

Code

代码语言:javascript
复制
#pragma GCC optimize("Ofast")
#pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,avx2,fma")
#pragma GCC optimize("unroll-loops")
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define W while
#define I inline
#define RI register int
#define LL long long
#define Cn const
#define CI Cn int&
using namespace std;
namespace Debug{
    Tp I void _debug(Cn char* f,Ty t){cerr<<f<<'='<<t<<endl;}
    Ts I void _debug(Cn char* f,Ty x,Ar... y){W(*f!=',') cerr<<*f++;cerr<<'='<<x<<",";_debug(f+1,y...);}
    Tp ostream& operator<<(ostream& os,Cn vector<Ty>& V){os<<"[";for(Cn auto& vv:V) os<<vv<<",";os<<"]";return os;}
    #define gdb(...) _debug(#__VA_ARGS__,__VA_ARGS__)
}using namespace Debug;
namespace FastIO{
    #define FS 100000
    #define tc() (FA==FB&&(FB=(FA=FI)+fread(FI,1,FS,stdin),FA==FB)?EOF:*FA++)
    #define pc(c) (FC==FE&&(clear(),0),*FC++=c)
    int OT;char oc,FI[FS],FO[FS],OS[FS],*FA=FI,*FB=FI,*FC=FO,*FE=FO+FS;
    I void clear() {fwrite(FO,1,FC-FO,stdout),FC=FO;}
    Tp I void read(Ty& x) {x=0;RI f=1;W(!isdigit(oc=tc())) f=oc^'-'?1:-1;W(x=(x<<3)+(x<<1)+(oc&15),isdigit(oc=tc()));x*=f;}
    Ts I void read(Ty& x,Ar&... y) {read(x),read(y...);}
    Tp I void write(Ty x) {x<0&&(pc('-'),x=-x);W(OS[++OT]=x%10+48,x/=10);W(OT) pc(OS[OT--]);}
    Tp I void writeln(Ty x) {x<0&&(pc('-'),x=-x);W(OS[++OT]=x%10+48,x/=10);W(OT) pc(OS[OT--]);pc('\n');}
}using namespace FastIO;
Cn int N=5e3+10,XX=1e9+7;
int n,k,a[N],sa[N],sb[N],ic[N],Ans;
I int QP(RI a,RI b,CI p=XX){RI s=1;W(b) b&1&&(s=1LL*s*a%p),a=1LL*a*a%p,b>>=1;return s;}
// class Poly{
//     private:
//         int P,L,R[N<<2],B[N<<2];
//         I void NTT(int *s,CI op){
//             RI i,j,k,x,y,U,S;for(i=0;i<P;i++) i<R[i]&&(swap(s[i],s[R[i]]),0);
//             for(i=1;i<P;i<<=1) for(U=QP(QP(3,op,X),(X-1)/(i<<1),X),j=0;j<P;j+=i<<1)
//             for(S=1,k=0;k<i;k++,S=1LL*S*U%X) s[j+k]=((x=s[j+k])+(y=1LL*S*s[i+j+k]%X))%X,s[i+j+k]=(x-y+X)%X;
//         }
//     public:
//         int X,A[N<<2];
//         I void Mul(CI n,int *a,CI m,int *b){
//             RI i,t;P=1,L=0;W(P<=n+m) P<<=1,++L;for(i=0;i<P;i++) R[i]=(R[i>>1]>>1)((i&1)<<L-1);
//             for(i=0;i<P;i++) A[i]=B[i]=0;for(i=0;i<=n;i++) A[i]=a[i];for(i=0;i<=m;i++) B[i]=b[i];
//             for(NTT(A,1),NTT(B,1),i=0;i<P;i++) A[i]=1LL*A[i]*B[i]%X;
//             for(t=QP(P,X-2,X),NTT(A,X-2),i=0;i<=n+m;i++) A[i]=1LL*A[i]*t%X; 
//         }
// }T[3];
// I LL CRT(LL r1,LL p1,LL r2,LL p2,CI fg){
//     LL k=1LL*((r2-r1)%p2+p2)*QP(p1%p2,p2-2,p2)%p2;
//     return fg?((p1%XX)*k+r1)%XX:(p1*k+r1)%(p1*p2);
// }
// I void Mul(CI n,int *a,CI m,int *b){
//     RI i;for(i=0;i<3;i++) T[i].Mul(n,a,m,b);
//     for(i=0;i<=n+m;i++) a[i]=CRT(CRT(T[0].A[i],T[0].X,T[1].A[i],T[1].X,0),1LL*T[0].X*T[1].X,T[2].A[i],T[2].X,1);
// }
int T[N];
I void Mul(CI n,int *a,CI m,int *b){
    RI i,j;for(i=0;i<=n+m;i++) T[i]=0;
    for(i=0;i<=n;i++) for(j=0;j<=m;j++) T[i+j]=(1LL*a[i]*b[j]%XX+T[i+j])%XX;
    for(i=0;i<=n+m;i++) a[i]=T[i];
}
int main(){
    freopen("calculate.in","r",stdin),freopen("calculate.out","w",stdout);
    RI i,t;for(read(n,k),Ans=i=1;i<=n;i++) read(a[i]),Ans=1LL*Ans*a[i]%XX;//T[0].X=998244353,T[1].X=469762049,T[2].X=1004535809;
    for(sa[0]=a[1],sa[1]=XX-1,i=2;i<=n;i++) sb[0]=a[i],sb[1]=XX-1,Mul(i-1,sa,1,sb);
    for(ic[n]=QP(QP(n,n),XX-2),i=n-1;~i;i--) ic[i]=1LL*ic[i+1]*n%XX;
    for(t=1,i=0;i<=n;i++) Ans=(XX-1LL*sa[i]*t%XX*ic[i]%XX+Ans)%XX,t=1LL*t*(k-i)%XX;
    return writeln(Ans),cerr<<clock()<<'\n',clear(),0;
}
本文参与 腾讯云自媒体同步曝光计划,分享自作者个人站点/博客。
原始发表:2022-02-10 ,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 作者个人站点/博客 前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • YbtOJ 976「母函数」随机减法
    • Solution
      • Code
      领券
      问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档