POJ-1322 Chocolate（生成函数）

跳转到题目

Chocolate

Time Limit: 2000MS		Memory Limit: 65536K
				Special Judge

Description

In 2100, ACM chocolate will be one of the favorite foods in the world. 

"Green, orange, brown, red...", colorful sugar-coated shell maybe is the most attractive feature of ACM chocolate. How many colors have you ever seen? Nowadays, it's said that the ACM chooses from a palette of twenty-four colors to paint their delicious candy bits. 

One day, Sandy played a game on a big package of ACM chocolates which contains five colors (green, orange, brown, red and yellow). Each time he took one chocolate from the package and placed it on the table. If there were two chocolates of the same color on the table, he ate both of them. He found a quite interesting thing that in most of the time there were always 2 or 3 chocolates on the table. 

Now, here comes the problem, if there are C colors of ACM chocolates in the package (colors are distributed evenly), after N chocolates are taken from the package, what's the probability that there is exactly M chocolates on the table? Would you please write a program to figure it out? 

Input

The input file for this problem contains several test cases, one per line. 

For each case, there are three non-negative integers: C (C <= 100), N and M (N, M <= 1000000). 

The input is terminated by a line containing a single zero. 

Output

The output should be one real number per line, shows the probability for each case, round to three decimal places.

Sample Input

5 100 2

0

Sample Output

0.625 

第一反应是概率DP，但是直接推复杂度是O(n*c)

又想到可以用矩阵快速幂优化，复杂度为O(c^3*logn)，一看数据范围发现比直接递推没有优化多少

看了一下网上这种做法的题解，发现都是卡时间过或者TLE

也发现由于精度要求低，所以当n大于某一常数时，概率只与奇偶有关

所以可以如此进行剪枝计算

但是本题也可以延伸：求取出n块巧克力后桌面上剩余m块巧克力时共有多少不同的取巧克力序列？

这样的话就必须使用生成函数（貌似以前好多题都是求方案数的）

解法在：国家集训队2009论文集-母函数的性质及应用 中

但是感觉需要深厚的数学功底才行，要不然会很容易算错

#include <cstdio>
#include <cstring>
#include <algorithm>

using namespace std;

const int MAXN=107;
const int MOD=100;

int c,n,m;
double tp1[MAXN],tp2[MAXN],tp1_[MAXN],tp2_[MAXN];//tp1[k]表示{[e^x-e^(-x)]/2}^m展开式中e^kx的系数，tp1_[k]表示其展开式中e^(-kx）的系数；tp2[k]、tp2_[k]同理
double p[MAXN],p_[MAXN];//p[k]表示{[e^x-e^(-x)]/2}^m*{[e^x+e^(-x)]/2}^(c-m)展开式中e^kx的系数，p_[k]则表示其展开式中e^(-kx)的系数

double quickPow(double a,int b) {
    double res=1;
    while(b>0) {
        if((b&1)==1) {
            res*=a;
        }
        a*=a;
        b>>=1;
    }
    return res;
}

double C(double res,int a,int b) {
    if(b>a-b) {
        b=a-b;
    }
    for(int i=1;i<=b;++i) {
        res*=(a-i+1.0)/i;
    }
    return res;
}

void getCoefficient() {//{[e^x-e^(-x)]/2}^m*{[e^x+e^(-x)]/2}^(c-m)展开式中e^kx和e^(-kx)的系数
    for(int k=0;k<=c;++k) {
        tp1[k]=tp1_[k]=tp2[k]=tp2_[k]=p[k]=p_[k]=0;
    }
    double sign,halfEm=quickPow(0.5,m);
    int exp;
    for(int i=0;i<=m;++i) {//计算前半部分表达式展开式的e^exp的系数
        sign=(i&1)==1?-1:1;
        exp=m-i-i;
        if(exp<0) {
            tp1_[-exp]+=sign*C(halfEm,m,i);
        }
        else {
            tp1[exp]+=sign*C(halfEm,m,i);
        }
    }
    halfEm=quickPow(0.5,c-m);
    for(int j=0;j<=c-m;++j) {//计算后半部分表达式展开式的e^exp的系数
        exp=c-m-j-j;
        if(exp<0) {
            tp2_[-exp]+=C(halfEm,c-m,j);
        }
        else {
            tp2[exp]+=C(halfEm,c-m,j);
        }
    }
    double fp,sp;//分别表示第一个式子和第二个式子的系数
    for(int i=-m;i<=m;++i) {//计算整个表达式展开式的e^exp的系数
        fp=i<0?tp1_[-i]:tp1[i];
        for(int j=-(c-m);j<=c-m;++j) {
            sp=j<0?tp2_[-j]:tp2[j];
            exp=i+j;
            if(exp<0) {
                p_[-exp]+=fp*sp;
            }
            else {
                p[exp]+=fp*sp;
            }
        }
    }
}

double solve() {
    getCoefficient();
    double ans=0,sign=(n&1)==1?-1:1;
    for(int k=1;k<=c;++k) {
        ans+=C(quickPow(1.0*k/c,n)*(p[k]+sign*p_[k]),c,m);
    }
    return ans;
}

int main() {
    while(scanf("%d",&c),c!=0) {
        scanf("%d%d",&n,&m);
        if(m>n||m>c||((n+m)&1)==1) {
            printf("0.000\n");
        }
        else if(n==0&&m==0){
            printf("1.000\n");
        }
        else {
            printf("%.3lf\n",solve());
        }
    }
    return 0;
}