原題是機器人走方格的問題：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。程式碼與上一種方法類似，在此省略。