[POJ 1322] Chocolate (生成函数 或 概率dp)

题目描述

Problem Description
In 2100, ACM chocolate will be one of the favorite foods in the world.

“Green, orange, brown, red…”, colorful sugar-coated shell maybe is the most attractive feature of ACM chocolate. How many colors have you ever seen? Nowadays, it’s said that the ACM chooses from a palette of twenty-four colors to paint their delicious candy bits.

One day, Sandy played a game on a big package of ACM chocolates which contains five colors (green, orange, brown, red and yellow). Each time he took one chocolate from the package and placed it on the table. If there were two chocolates of the same color on the table, he ate both of them. He found a quite interesting thing that in most of the time there were always 2 or 3 chocolates on the table.

Now, here comes the problem, if there are C colors of ACM chocolates in the package (colors are distributed evenly), after N chocolates are taken from the package, what’s the probability that there is exactly M chocolates on the table? Would you please write a program to figure it out?
Input
The input file for this problem contains several test cases, one per line.

For each case, there are three non-negative integers: C (C <= 100), N and M (N, M <= 1000000).

The input is terminated by a line containing a single zero.
Output
The output should be one real number per line, shows the probability for each case, round to three decimal places.

题目大意

一个口袋中装有巧克力，巧克力的颜色有$c$种。现从口袋中取出一个巧克力，若取出的巧克力与桌上某一已有巧克力颜色相同，则将两个巧克力都取走，否则将取出的巧克力放在桌上。设从口袋中取出每种颜色的巧克力的概率均等。求取出$N$个巧克力后桌面上剩余$M$个巧克力的概率。注意是无限输入，读入到$c=0$为止

样例输入输出

Sample Input
5 100 2
0

Sample Output
0.625

解题思路

生成函数

根据题意，只有当$m<=n$，$m<=c$ 且 $n-m$为偶数时才有解，下面我们只考虑有解的情况
首先我们可以根据题意将根据奇偶巧克力分成两类，其中，奇数个的有$m$种颜色，那么我们可以选择其中一种分类的方式，将其写成生成函数的形式，即：
$$(\frac{e^x-e^{-x}}{2}{)}^m\times (\frac{e^x+e^{-x}}{2}{)}^{c-m}$$
很容易知道最终的答案就为$\frac{a_n\times C_m^c\times n!}{c^n}$，其中$a_n$为上式展开后第n项的系数

---

最后的工作就是将上式展开，先二项式展开，得：
$$f(x)=2^{-c}\times (\sum _{i=0}^m(-1{)}^i\times C_m^i\times e^{(m-i)x}\times e^{-ix})\times (\sum _{j=0}^{c-m}{C}_{c-m}^j\times e^{jx}\times e^{-(c-m-j)x})$$
$$=2^{-c}\times \sum _{i=0}^m(-1{)}^i\times C_m^i\times e^{(m-2i)x}\times \sum _{j=0}^{c-m}{C}_{c-m}^j\times e^{(2j+m-c)x}$$
求一波卷积，得：
$$f(x)=2^{-c}\times \sum _{i=0}^m\sum _{j=0}^{c-m}(-1{)}^i\times C_m^i{C}_{c-m}^j\times e^{(2j-2i+2m-c)x}$$
再将$e^{(2i+2j+2m-c)x}$泰勒展开，得：
$$f(x)=2^{-c}\times \sum _{i=0}^m\sum _{j=0}^{c-m}(-1{)}^i\times C_m^i{C}_{c-m}^j\times \sum _{k=0}^{\infty }\frac{((2j-2i+2m-c)x{)}^k}{k!}$$

---

再重复一遍，最后答案为：
$$\frac{a_n\times C_m^c\times n!}{c^n}$$
代码里面没有$\times 2^{-c}$，故最后乘上；代码里把$\times c^{-n}$甩到了每一次$a_n$的累加里，故没有最后计算；代码里再求$a_n$时，没有算$\frac{1}{k!}$，故不用最后$\times n!$
（写的有一点点小混乱，主要是Markdown写公式源码很丑，虽然写好后长得很好看。。。）

概率dp

前面扯了一堆堆的高 qi 深 qi 莫 guai 测 guai 的玩意儿，相比之下，概率dp就简单易懂多了。。。
首先定义$dp[i][j]$表示桌子上已经摆出了$i$颗巧克力，最终剩下了$j$颗的概率，很容易给出转移方程式：
$$dp[i][j]=dp[i-1][j+1]\times \frac{j+1}{c}+dp[i-1][j-1]\times \frac{c-(j-1)}{c}$$
想要达到状态$(i,j)$有两种情况：

• 第一种，桌子上已有$j+1$种不同颜色的巧克力的各一颗，那么想要使其加上第$i$颗巧克力后只剩下$j$颗巧克力就需要拿出一颗与原本已有的巧克力之一的颜色相同的巧克力，然后吃掉，此事件发生的概率为$\frac{j+1}{c}$
• 第二种，桌子上只有$j-1$种不同颜色的巧克力，此时想要加上第$i$颗巧克力后剩下$j$颗就需要拿出一颗不同于桌上已有的巧克力的颜色的巧克力，此事件发生的概率为$\frac{c-(j+1)}{c}$

综合以上两种情况即可得到最上面给出的状态转移方程式

---

做法是给出来了，可是我们发现此算法的时间复杂度为$O(nm)$；回头看一眼数据范围：$C\le 100$ $N,M\le 1000000$，然后就可以 自闭 去了。。。
接着一波很迷的操作就来了，玄学得贼有道理。。。
我们发现数据范围很大，时间复杂度很大，但是，输出数据的精度很小，所以$N$越变越大，对答案的贡献越来越小，小到对答案产生不了影响，此时就可以停止计算了。这个时候，$N$大概在$500$上下（保险一点取到$1000$也不会$T$）

参考代码

生成函数

#include <cstdio>
#include <iostream>
using namespace std;

const int N = 100;
int c, n, m; double C[N + 5][N + 5], ans;

inline double qkpow (double x, int y) {
	double res = 1;
	for (; y; y >>= 1, x *= x)
		if (y & 1) res *= x; 
	return res;
}

inline void prepare () {
	C[0][0] = 1.0;
	for (int i = 1; i <= 100; ++ i) {
		C[i][0] = 1.0;
		for (int j = 1; j <= i; ++ j)
			C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
	}
}

int main () {
	prepare ();
	while (scanf ("%d", &c) && c) {
		scanf ("%d %d", &n, &m);
		if (m > c || m > n || (n - m) & 1) {
			printf ("0.000\n"); continue;
		}
		ans = 0.0;
		for (int i = 0; i <= m; ++ i)
			for (int j = 0; j <= c - m; ++ j)
				ans += (i & 1 ? -1 : 1) * C[m][i] * C[c - m][j] 
					   * qkpow ((2 * (m - i + j) - c) * 1.0 / c, n);
		ans = ans * C[c][m] / qkpow (2.0, c);
		printf ("%.3lf\n", ans);
	} 
	return 0;
}

概率dp
代码用的滚动数组，不滚动也可以

#include <cstdio>
#include <iostream>
#include <cstring>
using namespace std;

const int N = 100;
int c, n, m; double dp[2][N + 5];

int main () {
	while (scanf ("%d", &c) && c) {
		scanf ("%d %d", &n, &m);
		
		if (m > c || m > n || (n - m) & 1) {
			printf ("0.000\n"); continue;
		}
		if (n >= 500) n = 500 - (n & 1);
		
		memset (dp, 0, sizeof dp);
		dp[0][0] = dp[1][1] = 1.0;
		for (int i = 2; i <= n; ++ i) {
			dp[i & 1][0] = dp[(i - 1) & 1][1] / c;
			dp[i & 1][c] = dp[(i - 1) & 1][c - 1] / c;
			
			for (int j = 1; j <= i && j < c; ++ j)
				dp[i & 1][j] = (j + 1.0) / c * dp[(i - 1) & 1][j + 1]
							 + (c - j + 1.0) / c * dp[(i - 1) & 1][j - 1];
		}
		printf ("%.3f\n", dp[n & 1][m]);
		//不晓得为嘛,用 %lf 输出它会炸（听取WA声一片...）
	}
	return 0;
}