Chocolate(母函数)(POJ-1322)

题目

一个口袋中装有巧克力，巧克力的颜色有 $c$ 种。现从口袋中取出一个巧克力，若取出的
巧克力与桌上某一已有巧克力颜色相同，则将两个巧克力都取走，否则将取出的巧克力放在桌上。设从口袋中取出每种颜色的巧克力的概率均等。求取出 $n$ 个巧克力后桌面上剩余 $m$ 个巧克力的概率。
数据范围: $1\le c\le 100,n,m\le 1000000$

思路

我们首先考虑数据的合法性：
a. $n,m同奇偶$
b.$m<=nandm<=c$
那么我们可以发现 $n-m$ 为 $2$ 的倍数，那么剩下的 $m$ 颗意味着有 $m$ 种巧克力是取了奇数次，那我们就要确定到底是哪几种
然后我们考虑母函数表示
$$G(x)=C_c^m\ast (\frac{e^x-e^{-x}}{2}{)}^m\ast (\frac{e^x+e^{-x}}{2}{)}^{c-m}$$
解释一下,这表示我们首先确定取奇数的方案数 $C_c^m$ 然后由于取巧克力有顺序要求,所以我们考虑指数型母函数
于是我们要确定第 $n$ 项的系数
我们将 $e^x$ 保留进行还原多项式：
$$G(x)=2^{-c}\ast C_c^m\ast (e^x-e^{-x}{)}^m\ast (e^x+e^{-x}{)}^{c-m}$$
$$G(x)=2^{-c}\ast C_c^m\ast \sum _{i=0}^m(-1{)}^i\ast C_m^i\ast e^{(m-2i)x}\ast \sum _{j=0}^{c-m}{C}_{c-m}^j\ast e^{(c-m-2j)x}$$
$$G(x)=2^{-c}\ast C_c^m\ast \sum _{i=0}^m\sum _{j=0}^{c-m}(-1{)}^i\ast C_m^i{C}_{c-m}^j\ast e^{(c-2i-2j)x}$$
$$G(x)=2^{-c}\ast C_c^m\ast \sum _{i=0}^m\sum _{j=0}^{c-m}(-1{)}^i\ast C_m^i{C}_{c-m}^j\ast \sum _{k=0}^{+\infty }\frac{((c-2i-2j)x{)}^k}{k!}$$
那么第n项的系数 $a_n$ 为：
$$a_n=2^{-c}\ast C_c^m\ast \sum _{i=0}^m\sum _{j=0}^{c-m}(-1{)}^i\ast C_m^i{C}_{c-m}^j\ast \frac{(c-2i-2j{)}^n}{n!}$$
所以符合条件的方案数为 $a_n\ast n!$ ,我们发现总方案数为: $c^n$
所以概率为：
$$\frac{a_n\ast n!}{c^n}$$
注意,尽量将除法放在最后减小精度误差

代码

#include<set>
#include<map>
#include<stack>
#include<cmath>
#include<queue>
#include<deque>
#include<cstdio>
#include<ctime>
#include<bitset>
#include<vector>
#include<climits>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define LL long long
#define ULL unsigned long long
int read(){
	int f=1,x=0;char c=getchar();
	while(c<'0'||'9'<c){if(c=='-')f=-1;c=getchar();}
	while('0'<=c&&c<='9') x=(x<<3)+(x<<1)+c-'0',c=getchar();
	return f*x;
}
#define MAXN 100
#define INF 0x3f3f3f3f
double Pow(double x,int y){
	double ret=1.0;
	while(y){
		if(y&1) ret*=x;
		x*=x;
		y>>=1;
	}
	return ret;
}
double C[MAXN+5][MAXN+5];
void Prepare(){
	C[0][0]=1;
	for(int i=1;i<=MAXN;i++){
		C[i][0]=1;
		for(int j=1;j<=i;j++)
			C[i][j]=C[i-1][j-1]+C[i-1][j];
	}
	return ;
}
int main(){
	int c,n,m;
	Prepare();
	while(1){
		c=read();
		if(!c) break;
		n=read(),m=read();
		if(m>c||m>n||(n-m)%2){
			puts("0.000");
			continue;
		}
		double ans=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]*Pow((c-2*i-2*j)*1.0/c,n);
		ans=ans*C[c][m]/Pow(2,c);
		printf("%f\n",ans);
	}
	return 0;
}