【学习笔记】【UNR #3】百鸽笼

感觉自己对一些东西的理解还是不到位啊。

原问题可以转化成,有 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(ai1,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 t1次操作中恰好对 a 1 a_1 a1操作了 a 1 − 1 a_1-1 a11次,其他 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)=iti!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) (a11)!xa11i=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)=(a11)!xa11i=2n(yTai(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=(ni)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<<" "; } }

__EOF__

本文作者仰望星空的蚂蚁
本文链接https://www.cnblogs.com/cqbzly/p/17529952.html
关于博主:评论和私信会在第一时间回复。或者直接私信我。
版权声明:本博客所有文章除特别声明外,均采用 BY-NC-SA 许可协议。转载请注明出处!
声援博主:如果您觉得文章对您有帮助,可以点击文章右下角推荐一下。您的鼓励是博主的最大动力!
posted @   仰望星空的蚂蚁  阅读(35)  评论(0编辑  收藏  举报  
相关博文:
阅读排行:
· 全程不用写代码,我用AI程序员写了一个飞机大战
· MongoDB 8.0这个新功能碉堡了,比商业数据库还牛
· 记一次.NET内存居高不下排查解决与启示
· 白话解读 Dapr 1.15:你的「微服务管家」又秀新绝活了
· DeepSeek 开源周回顾「GitHub 热点速览」
历史上的今天:
2022-06-23 【学习笔记】拟阵
点击右上角即可分享
微信分享提示