求组合数的C++实现

#include <iostream>
using namespace std;

int com(int n,int r)
{
    int i,j,s=1;
    if(n-r<r) r=n-r;
    for(i=0,j=1;i<r;i++){
        s*=(n-i);
        while(j<=r && s%j==0){
            s/=j;
            j++;
        }
    }
    return s;
}

int main()
{
    int n,m;
    while(cin>>n>>m){
    cout<<com(n,m)<<endl;
    }
    return 0;
}

上面的代码只适合较小的n，经测试，当n=33 时，对于所有小于等于 n 的 m均能计算出值。n>33时，仅对于较小的m适用。

在网上找了找，学会了计算组合数的新方法

http://hi.baidu.com/sulipol/blog/item/9a83278990ba1b98a5c2723d.html

下面的代码则利用了 素数唯一分解定理 

N=P1^e1*P2^e2*...*Pr^er 其中Pi为素数因子，ei是正整数。

n是素数当且仅当 r=1 且 e1=1

组合数计算算法流程：1,素数打表　　2,幂整数拆分　　3,快速幂模

蓝色字为引用:

有一个小函数，求n！中含数x的几次方

int getx(int n,int x)//计算n!中质因子2的出现次数
{
    if(n==0)
        return 0;
    return n/x+get2(n/x);
}

int getx(int n,int x)//非递归写法
{
ans=0;
while(n)
{
   ans+=n/x;
   n/=x;
}
return ans;
}

算法解释：

60/5=12；12/5=2；所以可以肯定60！中含有5^14

第一次的60/5 我们找到了5^1倍数的所有数，

第二次12/5 相当于(60/5)/2==60/5^2,集计算60中含有5^2倍数的个数

…………以此类推

注意：由于n/(x^i-1)计算已经包含了n/x^i的n/(x^i-1)部分，所以n/x^i仅算n/x^i

个x，而非(n/x^i)*i个x，

2.拆分后对分式进行降幂运算

x1[i]-=(x2[i]+x3[i]) 可以证明x1[i]>=0

证明：由于C(n,m)结果必然为整数，即n!%((n-m)!*m!)==0

所以n！为(n-m)!*m!的整数倍，即pi^x1i%(pi^x2i*pi^x3i)==0

即pi^(x1[i]-x2[i]-x3[i])为整数

x1[i]-x2[i]-x3[i]>=0

证毕

#include <iostream>
#include <cstring>
#include <cmath>
#define SET(a) memset(a,0,sizeof(a))
#define M 150000
#define N 1000001
#define MOD 20110413
/* 
C(n,m)=n! / ((n-m!)*m!)
n!=2^x1*3^x2*5^x3***pi^xi 
*/

int p[M];
int x1[M];//n!
int x2[M];//(n-m)!
int x3[M];//m!

int pn; //numbers of the prime number
using namespace std;

int getprime(){
    int i,j,k=1;
    p[0]=2;
    for (i=3;i<M;i+=2){
        int flag=1;
        for(j=2;j<=sqrt(i);j++){
            if(i%j==0) {
                flag=0;
                break;
            }
        }
        if(flag==1) p[k++]=i;
    }
    return k;
}

int getx(int n,int p){
    int x=0;
    while (n){
        x+=n/p;
        n/=p;
    }
    return x;
}

int div(int n,int x[]){
    int i;
    for (i=0;i<pn && n>=p[i];i++){
        x[i]=getx(n,p[i]);
    }
    return i;
}

__int64 b_pow(int *arr,int n)
{
    int i;
    __int64 res=1,a,b,d;
    for(i=0;i<n;i++){
        if(arr[i]){
            a=p[i];
            b=arr[i];
            d=1;
            while(b){
                if(b%2) d=(a*d) % MOD;
                a=(a*a) % MOD;
                b>>=1;
            }
            res=(res*d) % MOD;
        }
    }
    return res;
}

__int64 C(int n,int m){
    int len,i;
    len=div(n,x1);
    div(n-m,x2);
    div(m,x3);
    for (i=0;i<M;i++)
        x1[i]-=(x2[i]+x3[i]);
    return b_pow(x1,len);
}

int main() {
    int n,m;
    __int64 ans;
    SET(p);
    pn=getprime();
    while(cin>>n>>m){    
        SET(x1);
        SET(x2);
        SET(x3); 
        ans=C(n,m);
        cout<<ans<<endl;
    }
    return 0;
}