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?(\frac{e^x?e^{?x}}{2}{)}^m?(\frac{e^x+e^{?x}}{2}{)}^{c?m}$$
解释一下,这表示我们首先确定取奇数的方案数 $C_c^m$ 然后由于取巧克力有顺序要求,所以我们考虑指数型母函数
于是我们要确定第 $n$ 项的系数
我们将 $e^x$ 保留进行还原多项式：
$$G(x)=2^{?c}?C_c^m?(e^x?e^{?x}{)}^m?(e^x+e^{?x}{)}^{c?m}$$
$$G(x)=2^{?c}?C_c^m?\sum _{i=0}^m(?1{)}^i?C_m^i?e^{(m?2i)x}?\sum _{j=0}^{c?m}{C}_{c?m}^j?e^{(c?m?2j)x}$$
$$G(x)=2^{?c}?C_c^m?\sum _{i=0}^m\sum _{j=0}^{c?m}(?1{)}^i?C_m^i{C}_{c?m}^j?e^{(c?2i?2j)x}$$
$$G(x)=2^{?c}?C_c^m?\sum _{i=0}^m\sum _{j=0}^{c?m}(?1{)}^i?C_m^i{C}_{c?m}^j?\sum _{k=0}^{+\infty }\frac{((c?2i?2j)x{)}^k}{k!}$$
那么第n项的系数 $a_n$ 为：
$$a_n=2^{?c}?C_c^m?\sum _{i=0}^m\sum _{j=0}^{c?m}(?1{)}^i?C_m^i{C}_{c?m}^j?\frac{(c?2i?2j{)}^n}{n!}$$
所以符合条件的方案数为 $a_n?n!$ ,我们发现总方案数为: $c^n$
所以概率为：
$$\frac{a_n?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;
}