POJ 1322 Chocolate(母函数)

题目链接：http://poj.org/problem?id=1322

题意：经过简单转化，我们需要解决一个这样的问题：每一次从C样东西里面随机选一种，一共选N次，求恰有M种东西被选了基数次的概率。

数据范围：C<=100, N,M<=10^6

思路一：简单DP

dp[m][n]表示选了n次之后恰好有m种东西出现奇数次的概率。

状态转移为dp[m][n+1]=dp[m+1][n]*(m+1)/C+dp[m-1][n]*(C-m+1)/C

复杂度为O(CN)，达到10^8级别，而且有多组数据，超时。

思路二：矩阵快速幂

将dp[0..C][n]写成一个C+1维列向量，记为v[n]，从思路一的转移方程我们可以得到v[n+1]和v[n]的关系：

v[n+1]=A*v[n]，其中A是一个(C+1)*(C+1)矩阵：

|     0       1/C      0       ...       ...     0  |

|  C/C        0     2/C      ...       ...     0  |

|     0   (C-1)/C   0     3/C       ...     0  |

|    ..................................................  |

|    0         ...      ...       ...        0  C/C |

|    0        ...       ...      ...       1/C   0  |

采用矩阵快速幂可以计算出v[N]=A^N*v[0]，复杂度为O(C^3*logN)，仍然达到了10^6级别，超时。

虽然这个方法是不对的，不过自己第一次写快速幂，还是记录一下代码。

/*
	PROG: POJ1322
	PROB:
*/

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

#define DEBUG 1
#define LOG(...) do { if (DEBUG) fprintf(stderr, __VA_ARGS__); } while(0)

#define MAXN 105
struct Matrix {
	double a[MAXN][MAXN];
	int N;
	Matrix(int N) : N(N) {
		memset(a, 0, sizeof(a));
	}
	Matrix operator*(const Matrix& rhs) const {
		Matrix ret(N);
		for (int i = 0; i < N; ++i)
		for (int j = 0; j < N; ++j)
		for (int k = 0; k < N; ++k)
			ret.a[i][j] += a[i][k]*rhs.a[k][j];
		return ret;
	}
};

Matrix Pow(const Matrix& A, int n) {
	if (n==0) {
		Matrix B(A.N);
		for (int i = 0; i < A.N; ++i)
			B.a[i][i] = 1.0;
		return B;
	}
	else if (n==1) return A;
	else if (!(n&1)) return Pow(A*A, n>>1);
	else return A*Pow(A*A, n>>1);
}

int main(void) {
	int C, N, M;
	while (scanf("%d", &C), C) {
		scanf("%d%d", &N, &M);
		if (M>C) {
			printf("0.000\n");
			continue;
		}
		Matrix A(C+1);
		double c = C;
		for (int i = 1; i <= C; ++i)
			A.a[i-1][i] = A.a[C+1-i][C-i] = i/c;
		Matrix B(C+1);
		B = Pow(A, N);
		printf("%.3f\n", B.a[M][0]);
	}
	return 0;
}

思路三：推公式

既然一般的递推方法都太慢了，那么只能推出公式了。

我们首先考虑一个这样的问题：每次从C种物品中选一个，一共选N次，使得第i种物品有a_i件，有多少种选法？

这个问题与原问题相比，少了一个限制条件。这个问题是很容易的：有N!/(a_1! a_2! ... a_C!)种选法。

一般来说这个公式可以直接用排列组合算出来，但是现在我们需要对a_i加一些附加条件，然后筛选出符合这些条件的结果，所以我们可以从另一个角度来看这个简化的问题。

从形式上看，这和指数母函数乘法是类似的：

若f(x)=e^x=1+x^1/1!+x^2/2!+x^3/3!+....，把(f(x))^C写成b_i*x^i/i!求和的形式，则

b_N为所有N!/(a_1! a_2! ... a_C!)之和，其中a_j之和为N。

如果我们希望选出a_1....a_M为奇数的那些，则可以在响应的f(x)的项中改变一下系数：对于前M个f(x)，偶数次幂的系数变成0；其余的f(x)，奇数次幂的洗漱变成0，这项乘出来之后就可以把满足条件的N!/(a_1! a_2! ... a_C!)筛选出来了。

到了这一步之后就基本可以做出来了。剩下的是写出具体的公式。

计算的时候要注意有一些地方可能会超过double范围，所以公式的具体形式要稍微变一下。

直接上代码：

/*
	PROG: POJ1322
	PROB: generating function
*/

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

#define DEBUG 1
#define LOG(...) do { if (DEBUG) fprintf(stderr, __VA_ARGS__); } while(0)

#define MAXN 105
double A[MAXN], B[MAXN];
int main(void) {
	int C, N, M;
	while (scanf("%d", &C), C) {
		scanf("%d%d", &N, &M);
		if (M>C||M>N||((N-M)&1)) {
			printf("0.000\n");
			continue;
		}

		memset(B, 0, sizeof(B));
		B[0] = 1.0;
		for (int i = 0; i < C; ++i) {
			memcpy(A, B, sizeof(A));
			memset(B, 0, sizeof(B));
			for (int j = 0; j <= i; ++j)
				if (i < M) B[j] -= A[j], B[j+1] += A[j];
				else B[j] += A[j], B[j+1] += A[j];
		}

		double ans = 0.0;
		double c = C, m = M;
		for (int i = 0; i <= C; ++i)
			ans += B[i]*pow((2*i-c)/c, N);
		
		for (int i = 1; i <= C-M; ++i)
			ans *= (m+i)/i;
		ans /= pow(2.0, C);
		printf("%.3f\n", ans);
	}
	return 0;
}