组合数C(n,m)的计算

n个互不相同的数的全排列是n！个。

一个有n个元素的集合的含有m个元素子集的个数为C(n,m)。

C(n,m)的计算方式：

1.公式：C(n,m) = n!/((n-m)! * m!)，在算法上较难实现，阶乘很快会爆long long

2.递推：C(n,m) = C(n-1,m-1) + C(n-1,m)

在算法上当然会采用第二种方式计算，而且因为C(n,m)本身值很大，所以大多数碰见它的情况会取模。

long long C(int i, int j){
	if(i == j) return c[i][j] = 1;
	if(j == 0) return c[i][j] = 1;
	if(j == 1) return c[i][j] = i;
	if(c[i][j]) return c[i][j];
	return c[i][j] = (C(i-1, j-1)%M+C(i-1, j)%M)%M;
}

详细解释及证明可见 百度百科：http://baike.baidu.com/view/2104972.htm

但是仅仅这样还是不够的，还有更优的方法，设我们要求的是C(n,m)%P

设欧拉函数为phi(i),表示1~i-1中与i互质的数的个数；

设pri[i]表示P的第i个素因子；

设tim(x,i)表示将i素因数分解后x的次数(x为素数)；

设ni(x)表示x在剩余系P中的逆元，ni(x) = x^(phi(x)-1)；

设fac(x) = x!/(a1*a2*……*an) a∈{x|P%x == 0}；

设res = fac(n)*ni(fac(m))*ni(fac(n-m))；

然后对于每个P的素因子pri[i]，做：

res=res*pow(pri[i],tim[x+y][i]-tim[x][i]-tim[y][i])%P；

最后C(n,m) = res。

这个方法的详细推导需要很多数学知识，有很多我也不太理解，所以就先单单写下模板了。

#include <cstdio>
#include <cstring>
#include <algorithm>
#define M 1000005
typedef long long L;
using namespace std;

int n, m, prm;
L P, phi, pri[101];
L fac[M], tim[M][101];
bool prime[M];

L pow(int x, int y){
    L t = 1, xx = x;
    while(y){
	if(y&1) t *= xx, t %= P;
	xx *= xx, xx %= P, y >>= 1;
    }
    return t;
}

void get_prime(){
    for(int i = 2; i*2 <= P; i++) prime[i*2] = 1;
	for(int i = 3; i*i <= P; i++){
  	    if(prime[i]) continue;
	    for(int j = i*i; j <= P; j += i*2)
		prime[j] = 1;
    }
}

void get_pri(){
    for(int i = 2; i <= P; i++){
        if(prime[i]) continue;
	    if(P%i == 0) pri[++prm] = i;
    }
    /*              another method, don't need prime[M]
    int x = P;
    for(int i = 2; i*i <= x; i++){
        if(x%i == 0) pri[++prm] = i; 
        while(x%i == 0) x /= i;
    }
    if(x != 1) pri[++prm] = x;
    */
}

void get_phi(){
    phi = P;
    for(int i = 1; i <= prm; i++)
        phi = phi/pri[i]*(pri[i]-1);
}

void get_fac_tim(){
    fac[0] = 1;
    for(int i = 1; i <= M; i++){
        L x = i;
        for(int j = 1; j <= prm; j++){
            tim[i][j] = tim[i-1][j];
            while(x%pri[j] == 0){
                x /= pri[j];
                tim[i][j]++;
            }
        }
        fac[i] = fac[i-1] * x % P;
    }
}

L C(int n, int m){
    L res = fac[n]*pow(fac[m],phi-1)%P*pow(fac[n-m],phi-1)%P;
    for(int i = 1; i <= prm; i++)
        res = res*pow(pri[i],tim[n][i]-tim[m][i]-tim[n-m][i])%P;
    return res;
}

int main()
{
	scanf("%d", &P);
	get_prime();
	get_pri();
	get_phi();
	get_fac_tim();
	while(scanf("%d %d", &n, &m) && n) printf("%lld\n", C(n, m));
	return 0;
}