组合数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
*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
[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;
}