• 【UOJ449】[集训队作业2018]喂鸽子


    【UOJ449】[集训队作业2018]喂鸽子

    题面

    UOJ

    题解

    考虑( ext{min-max})容斥,那么答案为

    [sum_{i=1}^n(-1)^{i+1}{nchoose i}f_i ]

    其中(f_i)表示选出(i)只鸽子然后其中第一只鸽子被喂饱的期望时间。
    那么(f_i=sum_{j}p_j imes j),其中(p_j)表示上述过程发生在时刻(j)的概率,
    这样子不好统计,我们换一种统计方法:
    (f_i=sum_jp'_j)(p'_j)表示在(j)时刻选出的鸽子没有一只喂饱的概率,也就是上述过程发生时刻大于等于(j+1)的概率。
    再定义一个数组(g_{c,s})表示(c)只鸽子单独拿出来,时刻(s)没有喂饱的的概率。
    那么

    [egin{aligned} p'_i&=sum_{j=0}^{i}g_{c,j} {ichoose j}ig(frac {c}{n}ig)^jig(frac {n-c}{n}ig)^{i-j} end{aligned} ]

    也就是说

    [egin{aligned} f_c&=sum_{igeq 1}p'_i\ &=sum_{igeq 1}sum_{j=0}^ig_{c,j} {ichoose j}ig(frac {c}{n}ig)^jig(frac {n-c}{n}ig)^{i-j}\ &=sum_{j=0}^{c(k-1)}g_{c,j}ig(frac cnig)^jsum_{igeq 0}{i+jchoose i}ig(frac{n-c}{n}ig)^i end{aligned} ]

    注意到后面那一坨如果令(x=frac {n-c}c)的话,就是(sum_{igeq 0}{i+jchoose j}x^i=(frac {1}{1-x})^{j+1}=(frac nc)^{j+1})
    那么(f_c=frac ncsum_{j=0}^{c(k-1)} g_{c,j})
    如果暴力的话可以考虑加入一个鸽子然后dp:
    (g_{c,s}=sum_{i}g_{c-1,s-i}(frac 1c)^i(frac {c-1}{c})^{s-i}{schoose i})
    NTT 优化做到总复杂度(O(n^2k(log n+log k)))

    考虑怎么把这个(log)去掉,令 EGF (g_c(x)=({sum_{i=0}^{k} frac {x^i}{i!}})^c)(实际上是(k-1)这里为了方便所以写成(k),下同),那么

    [g(x)=1+x+frac {x^2}{2!}+...+frac {x^k}{k!}\ g'(x)=g(x)-frac {x^k}{k!}\ egin{aligned} g'_c(x)&=cg_{c-1}(x)ig(g(x)-frac {x^k}{k!}ig)\ &=cig(g_{c}(x)-frac {x^k}{k!}g_{c-1}(x)ig) end{aligned}\ g_{c,s+1}=frac{c}{s+1}(g_{c,s}-frac {1}{k!}g_{c-1,s-k}) ]

    于是可以递推了,最后(g_{c,s})乘上一个(frac {1}{c^s})就是概率了。

    代码

    #include <iostream> 
    #include <cstdio> 
    #include <cstdlib> 
    #include <cstring> 
    #include <cmath> 
    #include <algorithm> 
    using namespace std; 
    const int Mod = 998244353; 
    int fpow(int x, int y) { 
    	int res = 1; 
    	while (y) {
    		if (y & 1) res = 1ll * res * x % Mod; 
    		x = 1ll * x * x % Mod; 
    		y >>= 1; 
    	} 
    	return res; 
    } 
    const int MAX_N = 55, MAX_M = MAX_N * 1005; 
    int g[MAX_N][MAX_M], ifc[MAX_M], fac[MAX_M]; 
    int f[MAX_N]; 
    int N, K; 
    int C(int n, int m) { return 1ll * fac[n] * ifc[m] % Mod * ifc[n - m] % Mod; } 
    int main () { 
    #ifndef ONLINE_JUDGE 
        freopen("cpp.in", "r", stdin); 
    #endif 
    	cin >> N >> K; 
    	for (int i = fac[0] = 1; i <= N * K; i++) fac[i] = 1ll * fac[i - 1] * i % Mod; 
    	ifc[N * K] = fpow(fac[N * K], Mod - 2); 
    	for (int i = N * K - 1; ~i; i--) ifc[i] = 1ll * ifc[i + 1] * (i + 1) % Mod; 
    	for (int i = 0; i < K; i++) g[1][i] = ifc[i]; 
    	for (int i = 2; i <= N; i++) {
    		g[i][0] = 1; 
    		for (int j = 1; j <= i * (K - 1); j++) { 
    			g[i][j] = g[i][j - 1]; 
    			if (j - K >= 0) g[i][j] = (g[i][j] - 1ll * ifc[K - 1] * g[i - 1][j - K]) % Mod; 
    			g[i][j] = (g[i][j] + Mod) % Mod; 
    			g[i][j] = 1ll * i * ifc[j] % Mod * fac[j - 1] % Mod * g[i][j] % Mod; 
    		} 
    	} 
    	
    	for (int i = 1; i <= N; i++) {
    		int t = 1ll * ifc[i] * fac[i - 1] % Mod; 
    		for (int pw = 1, j = 0; j <= i * (K - 1); j++, pw = 1ll * pw * t % Mod) 
    			g[i][j] = 1ll * g[i][j] * fac[j] % Mod * pw % Mod; 
    	} 
    	for (int i = 1; i <= N; i++) { 
    		for (int j = 0; j <= i * (K - 1); j++) f[i] = (f[i] + g[i][j]) % Mod; 
    		f[i] = 1ll * f[i] * N % Mod * ifc[i] % Mod * fac[i - 1] % Mod; 
    	} 
    	int ans = 0; 
    	for (int i = 1; i <= N; i++) { 
    		if (i & 1) ans = (ans + 1ll * C(N, i) * f[i]) % Mod; 
    		else ans = (ans - 1ll * C(N, i) * f[i]) % Mod, ans = (ans + Mod) % Mod; 
    	} 
    	printf("%d
    ", ans); 
        return 0; 
    } 
    
  • 相关阅读:
    GitHub里的Hello World!
    4 款消息队列软件产品大比拼(转)
    .net常用组件
    Dapper.NET使用(转)
    设置MYSQL允许用IP访问
    test1
    SQLServer 2008以上误操作数据库恢复方法——日志尾部备份(转)
    Quartz.NET配置
    Quartz CronTrigger配置
    Quartz CronTrigger最完整配置说明
  • 原文地址:https://www.cnblogs.com/heyujun/p/13039801.html
Copyright © 2020-2023  润新知