Choose and divide

Description
The binomial coefficient C(m,n) is defined as

                                   m!
             C(m,n) =  -----------
                               n!(m-n)!
Given four natural numbers p, q, r, and s, compute the the result of dividing C(p,q) by C(r,s).
The Input
Input consists of a sequence of lines. Each line contains four non-negative integer numbers giving values for p, q, r, and s, respectively, separated by a single space. All the numbers will be smaller than 10,000 with p>=q and r>=s.
The Output
For each line of input, print a single line containing a real number with 5 digits of precision in the fraction, giving the number as described above. You may assume the result is not greater than 100,000,000.
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
Output for Sample Input
0.12587
505606.46055
1.28223
0.48996
2.00000

3.99960

題意很簡單： 求C(p,q)/C(r,s);

思路：由於10^4，也不能for迴圈，去求C.所以先約分，利用唯一分解定理。

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<vector>
#define bug(a) cout<<a<<"--->\n";
using namespace std;

typedef unsigned long long ULL;

int prime[10001];
bool valid[10001];
int tot;
int e[10000];
int e_cnt;

typedef long long LL;

void get_prime()
{
    tot=0;
    memset(valid,true,sizeof valid);
    for(int i=2;i<10000;i++)
    {
        if(valid[i])
            prime[++tot]=i;
        for(int j=1;j<=tot&&i*prime[j]<10000;j++)
        {
            valid[i*prime[j]]=false;
            if(i%prime[j]==0) break;
        }
    }
}

void prime_break(int n,int d)
{

    for(int i=1;i<=tot;i++)
    {
        while(n%prime[i]==0)
        {
            n/=prime[i];
            e[i]+=d;
            if(i>e_cnt) e_cnt=i;
        }    
        if(n==1)
            return;
    }
}

void solve(int n,int d)
{
    for(int i=1;i<=n;i++)
        prime_break(i,d);
}

int main()
{
    get_prime();
    int p,q,r,s;
    while(scanf("%d%d%d%d",&p,&q,&r,&s)!=EOF)
    {
        memset(e,0,sizeof e);
        solve(p,1);
        solve(r-s,1);
        solve(s,1);
        solve(p-q,-1);
        solve(q,-1);
        solve(r,-1);
        double res=1.0;
        for(int i=1;i<=e_cnt;i++)
            if(e[i])
                res*=pow(prime[i],e[i]);

        printf("%.5lf\n",res);
    }

    return 0;
}