标签：mat 代码 验证 fct += clu 操作 scan 组合
老头子是小学校长，小学生（大哥）们都很听老头子的话。一天，老头子给小学生（大哥）们发苹果吃。
一共有 n 个小学生（大哥），老头子每一次会等概率选择一位小学生（大哥）并给他一个苹果。一个小学生（大哥）变得开心当且仅当他拥有的苹果数 ≥k。
因为老头子年纪大了，所以他想要你告诉他，期望多少次之后所有的小学生（大哥）都变得开心。
假设答案的最简分数形式为 a/b，你需要求出 w，满足 a≡b?w mod 998244353，且 0≤w<998244353。
Input format
一行两个数 n，k。
Output format
一行一个数表示答案。
Sample input
1 1
Sample output
1
Constraints
对于所有的数据，n≤50，k≤1000。
uoj#449. 【集训队作业2018】喂鸽子
你可以通过搜索上面那道题的题解完成此题。
嗯好我们进入正题。首先一眼 min-max 容斥。问题转化 n 个人中选，前 m 个人第一个拥有苹果数 ≥ k 的期望次数。
假如我们知道在 m 个人中选，第一个拥有苹果数 ≥ k 的期望次数为 f[m]。则 n 个人中选，前 m 个人第一个拥有苹果数 ≥ k 的期望次数为 f[m]*n/m。
为什么？因为期望等于概率的倒数！有恰当的正确性证明，只是我概率论学得不好。
考虑怎么求解 f[m] 。我们假设如果出现拥有苹果数 ≥ k 的事件后，强制结束操作。
设在发第 i 个苹果之前，操作还未结束的概率为 p[m][i]；设发了 i 个苹果，还是没有结束的方案数（即没有一个人的苹果数 ≥ k）为 g[m][i]；设发了 i 个苹果，刚好结束（即第 i 个苹果发给了之前某个拥有 k-1 个苹果的人，使那个人拥有了 k 个苹果）的方案数为 h[m][i]。显然以上定义中第二维的大小 <= m*k。
于是可以得到如下几个关系式：
\[p[m][i] = p[m][i-1]*\frac{g[m][i-1]}{g[m][i-1]+h[m][i-1]}\]
\[h[m][i] = m*C_{i-1}^{k-1}*g[m-1][i-k]\]
\[g[m][i] = \sum_{j=0}^{k-1}g[m-1][i-j]*C_{i}^{j}\]
\[f[m] = \sum_{i=0}^{m*k}\frac{h[m][i]}{g[m][i]+h[m][i]}*i\]
注意到 g 的转移式中，可以通过将组合数拆成阶乘形式来进行卷积。
于是时间复杂度为 O(n^2klog(nk))。
最终答案 ans 的表达式为：
\[ans = \sum_{i=1}^{n}f[i]*\frac{n}{i}*C_n^i*(-1)^{i-1}\]
代码实现的部分细节可能与上面描述的不一样（比如我的 p[m][i] 与 h[m][i] 是边算边求的），自行理解一下吧（因为我懒得改代码了）。
#include<cstdio>
#include<algorithm>
using namespace std;
const int G = 3;
const int MAXN = 50;
const int MAXK = 1000;
const int MOD = 998244353;
int pow_mod(int b, int p) {
int ret = 1;
while( p ) {
if( p & 1 ) ret = 1LL*ret*b%MOD;
b = 1LL*b*b%MOD;
p >>= 1;
}
return ret;
}
int pw[20 + 5], ipw[20 + 5];
int fct[MAXN*MAXK + 5], inv[MAXN*MAXK + 5];
void init() {
fct[0] = 1;
for(int i=1;i<=MAXN*MAXK;i++)
fct[i] = 1LL*fct[i-1]*i%MOD;
inv[MAXN*MAXK] = pow_mod(fct[MAXN*MAXK], MOD-2);
for(int i=MAXN*MAXK-1;i>=0;i--)
inv[i] = 1LL*inv[i+1]*(i + 1)%MOD;
for(int i=1;i<=20;i++)
pw[i] = pow_mod(G, (MOD-1)/(1<<i)), ipw[i] = pow_mod(G, (MOD-1) - (MOD-1)/(1<<i));
}
int comb(int n, int m) {
return 1LL*fct[n]*inv[m]%MOD*inv[n-m]%MOD;
}
void ntt(int *A, int len, int type) {
for(int i=0,j=0;i<len;i++) {
if( i < j ) swap(A[i], A[j]);
for(int l=(len>>1);(j^=l)<l;l>>=1);
}
for(int i=1;(1<<i)<=len;i++) {
int s = (1<<i), t = (s>>1);
int u = (type == 1) ? pw[i] : ipw[i];
for(int j=0;j<len;j+=s) {
for(int k=0,p=1;k<t;k++,p=1LL*p*u%MOD) {
int x = A[j+k], y = 1LL*A[j+k+t]*p%MOD;
A[j+k] = (x + y)%MOD, A[j+k+t] = (x + MOD - y)%MOD;
}
}
}
if( type == -1 ) {
int inv = pow_mod(len, MOD-2);
for(int i=0;i<len;i++)
A[i] = 1LL*A[i]*inv%MOD;
}
}
int dp[MAXN + 5][MAXN*MAXK + 5], A[2*MAXN*MAXK + 5], B[2*MAXN*MAXK + 5], C[2*MAXN*MAXK + 5];
int main() {
init();
int n, k, ans = 0;
scanf("%d%d", &n, &k);
dp[0][0] = 1;
for(int i=1;i<=n;i++) {
int len; for(len = 1; len <= i*(k-1); len <<= 1);
for(int j=0;j<len;j++)
A[j] = B[j] = 0;
for(int j=0;j<k;j++)
B[j] = inv[j];
for(int j=0;j<=(i-1)*(k-1);j++)
A[j] = 1LL*dp[i-1][j]*inv[j]%MOD;
ntt(A, len, 1), ntt(B, len, 1);
for(int j=0;j<len;j++)
C[j] = 1LL*A[j]*B[j]%MOD;
ntt(C, len, -1);
for(int j=0;j<=i*(k-1);j++)
dp[i][j] = 1LL*C[j]*fct[j]%MOD;
}
for(int i=1;i<=n;i++) {
int g = 0, p = 1;
for(int j=0;j<=(i-1)*(k-1);j++) {
int a = 1LL*i*dp[i-1][j]%MOD*comb(j+k-1, k-1)%MOD;
int b = dp[i][j+k];
int c = pow_mod(a+b, MOD-2);
g = (g + 1LL*p*a%MOD*c%MOD*(j+k)%MOD)%MOD;
p = 1LL*p*b%MOD*c%MOD;
}
int f = 1LL*pow_mod(i, MOD-2)*n%MOD*comb(n, i)%MOD*pow_mod(MOD-1, i-1)%MOD;
ans = (ans + 1LL*f*g%MOD)%MOD;
}
printf("%d\n", ans);
}这个题与 uoj#449. 【集训队作业2018】喂鸽子 可能的确是同一道题（连用暴力模拟求期望帮助验证正确性都是一样的）。
只是我不知道哪一道题先出现，所以也不好判断这个是否可以定义为原题。
好像是在雅礼集训的时候讲课讲过这个题，还讲过不用 min-max 容斥的做法。但是好像找不着当初那个课件。。。
考试后，我对比了一下我考场上想的做法和标算的做法，发现好像不太一样。。。
可以去搜 uoj 的那道题的题解了解一下正确的切题姿势。
注意刚好结束时，最后一颗苹果一定发给那个 k 个苹果的人。
就是因为这个细节一开始写挂了。。。
标签：mat 代码 验证 fct += clu 操作 scan 组合
原文地址：https://www.cnblogs.com/Tiw-Air-OAO/p/11146699.html
