uvaoj 10375 Choose and divide 唯一分解定理

uvaoj 10375 Choose and divide 唯一分解定理

已知C(m,n)=m!/(n!(m-n)!)，输入整数p,q,r,s（p>=q,r>=s,p,q,r,s<=10000)，计算C(p,q)/C(r,s)，保留5位小数。

可以使用唯一分解定理来做，用一个数组表示唯一分解式中各个素数的指数，可正可负。每乘上一个整数，就把相应的素因子指数上加上相应指数，除也一样。

代码如下：

/*************************************************************************
	> File Name: 10375.cpp
	> Author: gwq
	> Mail: gwq5210@qq.com 
	> Created Time: 2015年01月25日 星期日 10时34分28秒
 ************************************************************************/

#include <cmath>
#include <ctime>
#include <cctype>
#include <climits>
#include <cstdio>
#include <cstdlib>
#include <cstring>

#include <map>
#include <set>
#include <queue>
#include <stack>
#include <string>
#include <vector>
#include <sstream>
#include <iostream>
#include <algorithm>

#define INF (INT_MAX / 10)
#define clr(arr, val) memset(arr, val, sizeof(arr))
#define pb push_back
#define sz(a) ((int)(a).size())

using namespace std;
typedef set<int> si;
typedef vector<int> vi;
typedef map<int, int> mii;
typedef long long ll;
typedef unsigned long long ull;

const double esp = 1e-5;

#define N 10010

int primes[N], pcnt, isnpri[N], e[N];

void make_primes(void)
{
	pcnt = 0;
	clr(isnpri, 0);
	isnpri[0] = isnpri[1] = 1;
	for (int i = 2; i * i < N; ++i) {
		if (!isnpri[i]) {
			for (int j = i * i; j < N; j += i) {
				isnpri[j] = 1;
			}
		}
	}
	for (int i = 2; i < N; ++i) {
		if (!isnpri[i]) {
			primes[pcnt++] = i;
		}
	}
}

void add_int(int n, int d)
{
	for (int i = 0; i < pcnt; ++i) {
		while (n % primes[i] == 0) {
			n /= primes[i];
			e[i] += d;
		}
		if (n == 1) {		// 提前退出循环
			break;
		}
	}
}

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

int main(int argc, char *argv[])
{
	make_primes();
	int p, q, r, s;
	while (scanf("%d%d%d%d", &p, &q, &r, &s) != EOF) {
		clr(e, 0);
		add_fac(p, 1);
		add_fac(q, -1);
		add_fac(p - q, -1);
		add_fac(r, -1);
		add_fac(s, 1);
		add_fac(r - s, 1);
		double ans = 1.0;
		for (int i = 0; i < pcnt; ++i) {
			ans *= pow(primes[i], e[i]);
		}
		printf("%.5f\n", ans);
	}

	return 0;
}