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: POJ1322PROB:*/#include <cstdio>#include <cstring>using namespace std;#define DEBUG 1#define LOG(...) do { if (DEBUG) fprintf(stderr, __VA_ARGS__); } while(0)#define MAXN 105struct 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: POJ1322PROB: 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 105double 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;}