ZOJ1363 Chocolate

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

 

大致意思就是袋子里有C种巧克力，从里面取N块巧克力，如果有两个及以上颜色相同的就将它们全都吃掉，求桌上还剩几种颜色的巧克力

 

考虑递推式，就有P_i+1,k=P_i,k-1*(c-k+1)/c+P_i,k+1*(k+1)/c，然而对于这个方程来说，时间复杂度为Θ(nc)，跑不过

那不妨考虑矩阵将递推式转化为矩阵快速幂求解就将复杂度降下来了

构造矩阵Gⁿ*F，那么其中求得的解就是所求结果

具体实现如下

#include<bits/stdc++.h>
#define mem0(a) memset(a,0,sizeof(a))
using namespace std;
const int N=110;
struct Mat{
    int n;
    double m[N][N];
    Mat(int n){
        this->n=n;mem0(m);
    }
    Mat operator * (const Mat &b){
        Mat res(n);
        for(int i=0;i<n;++i)
            for(int j=0;j<n;++j)
                for(int k=0;k<n;++k)
                    res.m[i][j]+=m[i][k]*b.m[k][j];
        return res;
    }
};
Mat I(int n){
    Mat res(n);
    for(int i=0;i<n;++i)res.m[i][i]=1;
    return res;
}
Mat FPM(Mat a,int b){
    Mat res=I(a.n);
    for(;b;b>>=1,a=a*a)if(b&1)res=res*a;
    return res;
}
int c,n,m;
int main(){
    #ifndef ONLINE_JUDGE
        freopen("test.in","r",stdin);freopen("test.out","w",stdout);
    #endif
    while(~scanf("%d",&c)&&c){
        scanf("%d%d",&n,&m);
        if(!n&&!m){puts("1.000");continue;}
        if(m>n||m>c||(n+m)&1){puts("0.000");continue;}
        Mat ff(c+1);
        Mat gg(c+1);
        ff.m[0][0]=1;
        for(int i=0;i<=c;++i){
            if(i!=0)gg.m[i][i-1]=(double)(c-i+1)/c;
            if(i!=c)gg.m[i][i+1]=(double)(i+1)/c;
        }
        Mat ans=FPM(gg,n)*ff;
        printf("%.3lf\n",ans.m[m][0]);
    }
    return 0;
}