跳轉到題目

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塊巧克力時共有多少不同的取巧克力序列？

這樣的話就必須使用生成函式（貌似以前好多題都是求方案數的）

但是感覺需要深厚的數學功底才行，要不然會很容易算錯

#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;
}