题目链接: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;
}