感觉自己对一些东西的理解还是不到位啊。
原问题可以转化成,有 n n n个数 { a i } \{a_i\} {ai},每次从 n n n个数中随机选一个 i i i,令 a i = max ( a i − 1 , 0 ) a_i=\max(a_i-1,0) ai=max(ai−1,0),在期望有限步后 ∑ a i \sum a_i ∑ai会等于 1 1 1(可能超过 n n n步),设随机变量 X X X为此时唯一非零的位置的下标,求 X X X的分布。
再仔细思考一下。实际上这个转化和猎人杀是一个东西。就是说,对于不合法的操作就一直继续做直到合法,可以证明这不会影响最终的概率。具体怎么计算呢,考虑固定 a 1 a_1 a1为最后非零的位置,假设经过 t t t步后停下来了,那么对答案的贡献为 1 n t \frac{1}{n^t} nt1。注意这个 t t t是一个无穷级数。但是没关系,我们后面会说到怎么来处理它。
首先可以把答案用生成函数表达出来。有了前面那个转化的铺垫,我们可以任意的延长操作序列,这就很 n b nb nb了,我们可以延长成把 ∑ a i \sum a_i ∑ai变成 0 0 0时停止,但是要求最后一次操作的是 a 1 a_1 a1,那么就相当于在前 t − 1 t-1 t−1次操作中恰好对 a 1 a_1 a1操作了 a 1 − 1 a_1-1 a1−1次,其他 a i a_i ai至少操作了 a i a_i ai次。于是设 G t ( x ) = ∑ i ≥ t x i i ! G_t(x)=\sum_{i\ge t}\frac{x^i}{i!} Gt(x)=∑i≥ti!xi,那么答案的生成函数就是 x a 1 − 1 ( a 1 − 1 ) ! ∏ i = 2 n G a i ( x ) \frac{x^{a_1-1}}{(a_1-1)!}\prod_{i=2}^nG_{a_i}(x) (a1−1)!xa1−1∏i=2nGai(x)。注意可以用到 e x e^x ex的生成函数来替换,将 e x e^x ex看作变量 y y y就能得到一个 二元生成函数 ,因此设 T t ( x ) = ∑ i < t x i i ! T_{t}(x)=\sum_{i<t}\frac{x^i}{i!} Tt(x)=∑i<ti!xi,那么最终的生成函数 F ( x , y = e x ) = x a 1 − 1 ( a 1 − 1 ) ! ∏ i = 2 n ( y − T a i ( x ) ) F(x,y=e^x)=\frac{x^{a_1-1}}{(a_1-1)!}\prod_{i=2}^n(y-T_{a_i}(x)) F(x,y=ex)=(a1−1)!xa1−1∏i=2n(y−Tai(x))。这一步具体有什么用我们还暂时不清楚,不过似乎结构得到了简化(?)。
最后计算每一项 y i x j y^ix^j yixj对答案产生的贡献。将 e x i e^{xi} exi泰勒展开(终于用到之前写过的东西了),后面那个 x j x^j xj相当于就是平移,然后对应项乘上贡献的系数,写出来就是 j ! n j + 1 ∑ ( i n ) s ( s + j j ) \frac{j!}{n^{j+1}}\sum (\frac{i}{n})^s\binom{s+j}{j} nj+1j!∑(ni)s(js+j),注意看后面那个式子,我们 似乎见过 ,那么我就不写了,可以直接得到 L H S = j ! ( n − i ) j + 1 LHS=\frac{j!}{(n-i)^{j+1}} LHS=(n−i)j+1j!。这样收拢过后我们就可以解释换元的作用了,只需要关注 y i y^i yi前面的系数就好了,这也是生成函数的思想。
因为这题多项式次数比较小所以就暴力乘起来然后对每一项算就好了。算一项的复杂度大概是 O ( n 5 ) O(n^5) O(n5),如果每一项都暴力算的话就超时了。但是神奇的地方就在于 背包是可以撤销物品的 ,更进一步的说 多项式是可以做逆运算的,简单来说就是递推。感觉不太好用多项式求逆来解释所以就算了。
复杂度 O ( n 5 ) O(n^5) O(n5)。
#include<bits/stdc++.h>
#define ll long long
#define fi first
#define se second
#define pb push_back
#define db double
using namespace std;
const int mod=998244353;
const int N=35;
const int M=35*35;
int n,m,a[N],deg;
ll f[N][M],g[N][M],fac[M],inv[M],invnum[M];
void add(ll &x,ll y){x=(x+y)%mod;}
ll fpow(ll x,ll y=mod-2){
ll z(1);
for(;y;y>>=1){
if(y&1)z=z*x%mod;
x=x*x%mod;
}
return z;
}
void init(int n){
fac[0]=1;for(int i=1;i<=n;i++)fac[i]=fac[i-1]*i%mod;
inv[n]=fpow(fac[n]);for(int i=n;i>=1;i--)inv[i-1]=inv[i]*i%mod,invnum[i]=inv[i]*fac[i-1]%mod;
}
ll work(int i){
ll res(0);
memcpy(g,f,sizeof f);
for(int j=0;j<=n;j++){
for(int k=0;k<=deg;k++){
if(j)add(g[j][k],-g[j-1][k]);
for(int l=1;l<a[i];l++){
if(k>=l)add(g[j][k],g[j][k-l]*inv[l]);
}
g[j][k]=-g[j][k];
}
}
for(int j=0;j<n;j++){
ll mul=invnum[n-j];
for(int k=0;k<=deg;k++){
if(k>=a[i]-1)add(res,g[j][k-a[i]+1]*inv[a[i]-1]%mod*fac[k]%mod*mul);
mul=mul*invnum[n-j]%mod;
}
}
return res;
}
signed main(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
cin>>n;for(int i=1;i<=n;i++)cin>>a[i],deg+=a[i];
f[0][0]=1;init(deg);
for(int i=1;i<=n;i++){
memset(g,0,sizeof g);
for(int j=0;j<i;j++){
for(int k=0;k<=deg;k++){
if(f[j][k]){
add(g[j+1][k],f[j][k]);
for(int l=0;l<a[i];l++){
add(g[j][k+l],-f[j][k]*inv[l]);
}
}
}
}
memcpy(f,g,sizeof g);
}
for(int i=1;i<=n;i++){
cout<<(work(i)+mod)%mod<<" ";
}
}