POJ 2613 Choose and divide（数论）

Description
给出p,q,r,s，输出C(p,q)/C(r,s)
Input
多组用例，每组用例占一行包括四个整数p,q,r,s(p>=q,r>=s)，以文件尾结束输入
Output
对于每组用例，输出C(p,q)/C(r,s)，结果保留到小数点后五位
Sample Input
10 5 14 9
93 45 84 59
145 95 143 92
995 487 996 488
2000 1000 1999 999
9998 4999 9996 4998
Sample Output
0.12587
505606.46055
1.28223
0.48996
2.00000
3.99960
Solution
阶乘的质因数分解，n!中素因子p的幂次为[n/p]+[[n/p]/p]+…（[x]表示对x下取整）
Code

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>d;
#define maxn 11111
typedef long long ll;
int prime[maxn],is_prime[maxn],res,num[maxn];
void get_prime()
{
    memset(is_prime,0,sizeof(is_prime));
    res=0;
    for(int i=2;i<maxn;i++)
        if(!is_prime[i])
        {
            prime[res++]=i;
            for(int j=i;j<maxn;j+=i)
                is_prime[j]=1;
        }
} 
void divide(int n,int state)
{
    for(int i=0;i<res;i++)
    {
        int m=n;
        while(m)
        {
            num[i]+=m/prime[i]*state;
            m/=prime[i];
        }
    }
}
int main()
{
    get_prime();
    int p,q,r,s;
    while(~scanf("%d%d%d%d",&p,&q,&r,&s))
    {
        memset(num,0,sizeof(num));
        double ans=1;
        if(p>q)divide(p,1),divide(q,-1),divide(p-q,-1);
        if(r>s)divide(r,-1),divide(s,1),divide(r-s,1);
        for(int i=0;i<res;i++)
            if(num[i]>0)ans*=pow(1.0*prime[i],1.0*num[i]);
            else if(num[i]<0)ans/=pow(1.0*prime[i],-1.0*num[i]);
        printf("%.5lf\n",ans);
    }
    return 0;
}