[生成函数阶段性小结][CF891E]Lust

问题描述

给你一个长度为n的数组a[]，还有操作数K，每次操作你在下标[1..n]中等概率选择一个下标x，贡献+=∏i!=xa[i]∏i!=xa[i]\prod_{i!=x}a[i]，然后a[x]-=1。求K次操作后贡献期望值，对1e9+7取模。
n<=5000,K<=1e9

问题分析

首先期望可以看作是所有方案的贡献除以方案数。
我们看到那个∏∏\prod不能算x，这不利于我们化式子，直接这样根本没有思路。
我们稍微化一化：∏ia[i]−∏ia′[i]∏ia[i]−∏ia′[i]\prod_{i}a[i]-\prod_{i}a'[i]，其中a’为a[x]-=1后的序列。
我们观察连续几次相减，发现中间的项全部被消掉了。
那么设b[i]表示i这个位置被删了b[i]次。那么一种方案的贡献就是∏a[i]−∏(a[i]−b[i])∏a[i]−∏(a[i]−b[i])\prod{a[i]}-\prod{(a[i]-b[i])}。
那么E(贡献)=∏ai−K!∏(ai−bi)nkE(贡献)=∏ai−K!∏(ai−bi)nkE(贡献)=\prod{a_i}-\frac{K!\prod(a_i-b_i)}{n^k}。
其中b[i]要枚举所有情况。K!表示，当b[]确定时，每次操作的x的顺序有K!。
那么问题变成了怎么把所有b[]的情况的贡献计算出来。

生成函数

一般生成函数

记A(x)=∑n≥0aixiA(x)=∑n≥0aixiA(x)=\sum_{n≥0}{a_ix^i}为原序列a[]的一个生成函数，由于一般使用生成函数，我们是用来化简一些转移的，所以并不关心x的具体值是什么，我们只关注系数。
就比如背包问题，把记录当前做到的每种体积的方案的数组f[]变成生成函数，我们和一个物品的生成函数相乘，即可得到转移过后的f[]的生成函数。
上面这种转移我们从多项式的角度理解，即卷积，写过FFT的人应该都能理解。
生成函数的笛卡尔积就是两个生成函数A和B的卷积。
生成函数的优美在于，可以把一些特殊的形式幂级数化为一个低次的式子.
比如A(x)=1+x+x2+...x+∞=11−xA(x)=1+x+x2+...x+∞=11−xA(x)=1+x+x^2+...x^{+∞}=\frac{1}{1-x}，这是等比数列求和，由于正无穷的项的系数我们不会用到，所以忽略了。
这样，几个多项式相乘的形式就可能很简单，也是我们做题的关键，即把答案的生成函数做出来，然后就可以展开得出系数，这道题就是一个应用。

几个重要的式子

OGF的：
11−x=∑i≥0xi11−x=∑i≥0xi\frac{1}{1-x}=\sum_{i≥0}{x^i}
1(1−ax)m=∑i≥0Cm−1i+m−1aixi1(1−ax)m=∑i≥0Ci+m−1m−1aixi\frac{1}{(1-ax)^m}=\sum_{i≥0}{C_{i+m-1}^{m-1}a^ix^i}
这个怎么理解呢？考虑一个背包问题，组合数的意义是i的体积可以由m种球共i个（体积为1）组成；a为某个系数，其实可以把a和x放在一起，就跟第一个式子一样了。
特别地，x(1−x)2=∑i≥1i∗xix(1−x)2=∑i≥1i∗xi\frac{x}{(1-x)^2}=\sum_{i≥1}{i*x^i}
分子的x相当于把原序列右移一位。
EGF（指数型）的，也就是形如A(x)=∑i≥0aii!xiA(x)=∑i≥0aii!xiA(x)=\sum_{i≥0}\frac{a_i}{i!}x^i，实际操作的时候我们储存的系数是aii!aii!\frac{a_i}{i!}
ex=∑i≥0xii!ex=∑i≥0xii!e^x=\sum_{i≥0}\frac{x^i}{i!}
ln(1+x)=∑i>0(−1)i+1xiiln(1+x)=∑i>0(−1)i+1xiiln(1+x)=\sum_{i>0}\frac{(-1)^{i+1}x^i}{i}
−ln(1−x)=∑i>0xii−ln(1−x)=∑i>0xii-ln(1-x)=\sum_{i>0}\frac{x^i}{i}
记住就好。后面两个没有0可以理解为不能被0除。
十分脑洞地，我们上面的x可以代表一个多项式，那么就可以出现非常好玩的情况。
比如A(x)=∑i≥0xii!=exA(x)=∑i≥0xii!=exA(x)=\sum_{i≥0}\frac{x^i}{i!}=e^x，B(x)=∑i>0(−1)i+1xii=ln(1+x)B(x)=∑i>0(−1)i+1xii=ln(1+x)B(x)=\sum_{i>0}\frac{(-1)^{i+1}x^i}{i}=ln(1+x)，那么A(x)B(x)=1+xA(x)B(x)=1+xA(x)^{B(x)}=1+x。
EGF有一个模型是带标号的对象拼接。
考虑一个带标号对象B是由若干个带标号对象A组合而成，A的个数不限制，那么就有B(x)=∑i≥0A(x)ii!=eA(x)B(x)=∑i≥0A(x)ii!=eA(x)B(x)=\sum_{i≥0}\frac{A(x)^i}{i!}=e^{A(x)}，大小为n的对象B的方案仍然是[xn]B(x)[xn]B(x)[x^n]B(x)。
除以阶乘的原因是，A之间是不可区分的，每种i个A的拼接方案算重i!次。
有了这些知识，我们可以轻松解决这道题了！

回到题目

E(贡献)=∏ai−K!∏(ai−bi)nkE(贡献)=∏ai−K!∏(ai−bi)nkE(贡献)=\prod{a_i}-\frac{K!\prod(a_i-b_i)}{n^k}
我们只用算后面那一块。
首先算(ai-bi)那块。
我们先分析某个位置i，我们设生成函数Fi(x)Fi(x)F_i(x)代表a[i]减去某个b[i]的贡献，我们要表示他的全部情况，不妨设Fi(x)=∑j≥0ai−jj!xjFi(x)=∑j≥0ai−jj!xjF_i(x)=\sum_{j≥0}\frac{a_i-j}{j!}x^j，第j次的系数就是b[i]=j这一位的贡献。
推式子
Fi(x)=∑j≥0ai−jj!xjFi(x)=∑j≥0ai−jj!xjF_i(x)=\sum_{j≥0}\frac{a_i-j}{j!}x^j
=∑j≥0aij!xj−1(j−1)!xj=∑j≥0aij!xj−1(j−1)!xj=\sum_{j≥0}\frac{a_i}{j!}x^j-\frac{1}{(j-1)!}x^j
注意到ai是常数。
=(ai−x)ex=(ai−x)ex=(a_i-x)e^x
那么所有的贡献和就是所有生成函数乘起来。
Ftot=enx∏i(ai−x)Ftot=enx∏i(ai−x)F_{tot}=e^{nx}\prod_{i}(a_i-x)。
分别考虑左右某一项的系数，左边是熟知的形式，右边可以暴力卷积，（像背包一样做）。设右边的第i项系数为c[i]，那么答案变成∑i=0..nc[i]∗nK−ii!∑i=0..nc[i]∗nK−ii!\sum_{i=0..n} c[i]*\frac{n^{K-i}}{i!}，i到n为止就行了，c只有0~n有东西。
带回去，∑i=0..nc[i]∗Ki↓ni，Ki↓表示∏j=K..K−i+1j∑i=0..nc[i]∗Ki↓ni，Ki↓表示∏j=K..K−i+1j\sum_{i=0..n} c[i]*\frac{K^{i↓}}{n^i}，K^{i↓}表示\prod_{j=K..K-i+1}j。
那么答案就出来了。

代码

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<set>
using namespace std;
typedef long long ll;
typedef double db;
#define fo(i,j,k) for(i=j;i<=k;i++)
#define fd(i,j,k) for(i=j;i>=k;i--)
#define cmax(a,b) (a=(a>b)?a:b)
#define cmin(a,b) (a=(a<b)?a:b)
const int N=5e3+5,mo=1e9+7,M=2e5+5;
int read()
{int x=0;char ch=getchar();while (ch<'0'||ch>'9') ch=getchar();while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}return x;
}
int n,K,a[N],F[N],G[N],rev[N],fac[N],ans,i,j,prod,In;
int ksm(int x,int y)
{int ret=1;while (y){if (y&1) ret=1ll*ret*x%mo;y>>=1;x=1ll*x*x%mo;}return ret;
}
int dwm(int x,int y)
{int ret=1;while (y--){ret=1ll*ret*x%mo;x--;}return ret;
}
int main()
{freopen("t4.in","r",stdin);//freopen("t4.out","w",stdout);scanf("%d %d",&n,&K);fo(i,1,n) scanf("%d\n",a+i);F[0]=1;fo(i,1,n){fd(j,i,1)F[j]=(1ll*F[j]*a[i]-F[j-1])%mo;F[0]=1ll*F[0]*a[i]%mo;}In=1;fo(i,0,n){ans=(ans+1ll*F[i]*dwm(K,i)%mo*ksm(In,mo-2))%mo;In=1ll*In*n%mo;}prod=1;fo(i,1,n) prod=1ll*prod*a[i]%mo;ans=(prod-ans)%mo;printf("%d\n",(ans+mo)%mo);
}