求大数n,m下组合数C(n+m,m)%Mod

原题是机器人走方格的问题：M * N的方格，一个机器人从左上走到右下，只能向右或向下走。有多少种不同的走法？由于方法数量可能很大，只需要输出Mod 10^9 + 7的结果。
此问题很简单，就直接是C(M+N-2,M-1)即可，但是当M+N很大时，是无法直接求出C(M+N-2,M-1)的，所以专门总结了一下大数下组合数的求解方法。

求大数的阶乘，如果不精确的话可以用斯特林公式：

1、转换为对数式，lg(C(n+m,m))=lg((n+m)!/(n!*m!))=lg((n+m)!)-lg(n!)-lg(m!)=lg((m+1)(m+2)…(m+n))-lg(n!)，将乘积化为加法运算。

//  对数法求C(n+m,m)
#include <iostream>
#include <cmath>
using namespace std;
typedef long long ll;
#define Mod 1000000009

ll Cnm(int n, int m)
{
    double sum1=0, sum2=0;
    ll res=(ll)1;
    for(int i=1; i<=m; i++)
    {
        sum1+=(double)log(i);
        sum2+=(double)log(n-m+i);
    }
    return (ll)exp(sum2-sum1)%Mod;
}

int main()
{
    int n,m;
    cin >> n >> m;
    cout << Cnm(n+m,(m>n?n:m)) << endl;
    system("pause");
    return 0;
}

但是此方法依然不适用与n、m过大的情况，但比直接计算阶乘适用范围要广得多。
2、同样不适用n过大情况的大数组合计算方法还有利用杨辉三角公式求解，C(n+m,m)=C(n+m-1,m-1)+C(n+m-1,m)，此公式证明很简单，展开拆分一项即可。

//  杨辉三角公式求C(n+m,m)
#include <iostream>
#include <vector>
using namespace std;
#define Mod 1000000009
typedef long long ll;

ll Cnm(int n, int m)
{
    vector<vector<ll>> Combi(n+1,vector<ll>(m+1,(ll)0)); 
    for(int i=0; i<(int)Combi[0].size(); i++)
        Combi[0][i]=(ll)0;
    for(int i=0; i<(int)Combi.size(); i++)
        Combi[i][0]=(ll)1;
    for(int i=1; i<=n; i++)
    {
        for(int j=1; j<=(i>m?m:i); j++)
        {
            Combi[i][j]=(Combi[i-1][j-1]+Combi[i-1][j])%Mod;
        }
    }
    return Combi.back().back();
}

int main()
{
    int n,m;
    cin >> n >> m;
    cout << Cnm(n+m,m) << endl;
    system("pause");
    return 0;
}

3、利用n!=2^a1*3^a2*…p^ak(p为质数),将原方程(n+m)!/(n!*m!)化为2^(a1-b1-c1)*3^(a2-b2-c2)…*p^(ak-bk-ck)，n!分解后p的次数为：n/p+n/p^2+…+n/p^k，则原式可以化成这k项式子对Mod取余的乘积。

//  质数化简求C(n+m,m)
#include <iostream>
#include <vector>
using namespace std;
#define Mod 1000000009
typedef long long ll;
//  计算n以内所有的质数
vector<int> primelessthanN(int n)
{
    vector<bool> isprime(n+1, true);
    vector<int> prime;
    prime.push_back(2);
    int i;
    for(i=3; i*i<=n; i+=2)
    {
        if(isprime[i])
        {
            prime.push_back(i);
            for(int j=i*i; j<=n; j+=i)
                isprime[j]=false;
        }
    }
    while(i<=n)
    {
        if(isprime[i])
            prime.push_back(i);
        i+=2;
    }
    return prime;
}
//  计算素数因子p的指数
int Cal(int n, int p)
{
    int res=0;
    ll rec=p;
    while((ll)n>=rec)
    {
        res+=(int)((ll)n/rec);
        rec*=(ll)p;
    }
    return res;
}
//  计算n^k对Mod取余
ll PowerMod(int n, int k)
{
    int t(k),n0(n);
    ll res=(ll)1;
    while(t)
    {
        if(t&1)
            res=(res*(ll)n0)%Mod;
        n0=(n0*n0)%Mod;
        t>>=1;
    }
    return res;
}
//  计算C(n,m)，即原式的C(n+m,m)
ll Cnm(int n, int m)
{
    vector<int> prime=primelessthanN(n);
    ll res=1;
    for(int i=0; i<prime.size(); i++)
    {
        res=(res*(PowerMod(prime[i],Cal(n,prime[i])-Cal(n-m,prime[i])-Cal(m,prime[i]))))%Mod;
    }
    return res;
}
//  main函数
int main()
{
    int n,m;
    cin >> n >> m;
    cout << Cnm(n+m,(m>n?n:m)) << endl;
    system("pause");
    return 0;
}

4、Lucas定理，p是一个素数，将n、m均转化为p进制数，表示如下：
（适用范围：n和m比较大，p是素数且比较小的时候）
m=m0+m1*p+m2*p^2+…mk*p^k,n=n0+n1*p+n2*p^2+…+nk*p^k，则C(n,m)=C(n0,m0)*C(n1,m1)*C(n2,m2)…*C(nk,mk)%p，即Lucas(n,m,p)=C(n%p,m%p)*Lucas(n/p,m/p,p)。
证明方法是利用p-进制的唯一性，具体证明如下

//  Lucas定理求C(n+m,m)
#include <iostream>
#include <vector>
using namespace std;
#define Mod 10009
typedef long long ll;
vector<long long> ff(Mod+1,1); 
//  预处理阶乘数组Combi
void InitFF(int mod)
{
    for(int i=1; i<=mod; i++)
        ff[i]=(ff[i-1]*i)%mod;
}
//  求最大公因数
int gcd(int a, int b)
{
    if(b==0) return a;
    return gcd(b,a%b);
}
//  阶线性同余方程，利用欧几里得定理，即求组合数分母阶乘式的
//  乘法逆元，具体看代码。
void _gcd(int a, int b, ll& x, ll& y)
{
    if(b==0)
    {
        x=(ll)1;
        y=(ll)0;
    }
    else
    {
        _gcd(b,a%b,x,y);
        ll temp=x;
        x=y;
        y=temp-(a/b)*y; //  x，y分别为ax+by=1的解，此方程的x满足(a*x)%y=1,即x为a的乘法逆元
        //  其实由费马小定理可知a的乘法逆元是a^(y-2)
    }
}
//  计算小n（在Mod范围内）时的组合数
ll C(int n, int m)
{
    if(n<m) return 0;
    ll res=(ll)1, x, y;
    int b=(int)(ff[n-m]*ff[m])%Mod;
    int a=(int)ff[n];
    int c=gcd(a,b);
    a/=c;
    b/=c;   //  保证a、b互质
    _gcd(b, Mod, x, y);
    x=(x+Mod)%Mod;      //  防止x为负数
    res=(x*a)%Mod;      //  求出C(n,m)
    return res;
}
//  Lucas函数
ll Cnm(int n, int m)
{
    InitFF(Mod);
    ll res=(ll)1;
    int m0(m), n0(n);
    while(m0||n0)
    {
        res=(res*C(n0%Mod, m0%Mod))%Mod;
        n0/=Mod;
        m0/=Mod;
    }
    return res; 
}

int main()
{
    int n,m;
    cin >> n >> m;
    cout << Cnm(n+m,(m>n?n:m)) << endl;
    system("pause");
    return 0;
}

这里简单介绍一下扩展欧几里得算法，是利用ax+by=c，其中c%gcd(a,b)=0，来求解满足条件的x,y值，其中x就是a对模b的乘法逆元，而y就是b对模a的乘法逆元，详解请参考http://blog.csdn.net/qq_27599517/article/details/50888092。
5、利用乘法逆元求解
此方法其实与上一种相同，求解乘法逆元除了用扩展欧几里得算法外，还可以直接由费马小定理得到（若对于a，p存在x，使得a*x mod p=1，那么我们称x为a对p的乘法逆元）,而费马小定理是：假如p是质数，且Gcd(a,p)=1，那么 a(p-1)（mod p）≡1。即：假如a是整数，p是质数，且a,p互质(即两者只有一个公约数1)，那么a的(p-1)次方除以p的余数恒等于1。
证明：
假如inv是b对于p的乘法逆元，即b*inv=p*t+1（t为整数），移项得inv=(p*t+1)/b
(a*inv) mod p
=(a*((p*t+1)/b)) mod p
=(a*(p*t/b+1/b)) mod p
=(a/b) mod p+(a*p*t/b) mod p
∵ (a*p*t/b) mod p=0
∴ 原式=(a/b) mod p
即(a*inv) mod p=(a/b) mod p
那么对于C(n+m)=(n+m)!/(n!*m!)就可以化为(n+m)!(n!*m!)^(p-2) mod p。代码与上一种方法类似，在此省略。