1. 程式人生 > >[組合數]求組合數的幾種方法總結

[組合數]求組合數的幾種方法總結

求C(n,m)%mod的方法總結

1.當n,m都很小的時候可以利用楊輝三角直接求。
C(n,m)=C(n-1,m)+C(n-1,m-1);

2.利用乘法逆元。
乘法逆元:(a/b)%mod=a*(b^(mod-2)) mod為素數。
逆元可以利用擴充套件歐幾里德或尤拉函式求得:

 1).擴充套件歐幾里德:b*x+p*y=1 有解,x就是所求

 2).費馬小定理:b^(p-1)=1(mod p),故b*b^(p-2)=1(mod p),所以x=b^(p-2)

1. n!/(m!*(n-m)! =x%p ,先對算出n!、m!、(n-m)!對p取模的餘數,就轉換為a/b=x%p;因為p為素數,所以等價於bx+py=a;然後用擴充套件的歐幾里得定理算出 bx’+py’=1的解,x=x’*a,就得到了最終的x的值,即C(m,n)%p得值。

2.逆元 其實如果mod是素數 則b的逆元其實就是b^(mod-2)
也就是 (m!(n-m)!)的逆元為 (m!(n-m)!)p-2 ;

int inv(int a) {  
    //return fpow(a, MOD-2, MOD);  
    return a == 1 ? 1 : (long long)(MOD - MOD / a) * inv(MOD % a) % MOD;  
}  
LL C(LL n,LL m)  
{  
    if(m < 0)return 0;  
    if(n < m)return 0;  
    if(m > n-m) m = n-m;  

    LL up = 1
, down = 1; for(LL i = 0 ; i < m ; i ++){ up = up * (n-i) % MOD; down = down * (i+1) % MOD; } return up * inv(down) % MOD; }

3.當n和m比較大,mod是素數且比較小的時候(10^5左右),通過Lucas定理計算

Lucas定理:A、B是非負整數,p是質數。A B寫成p進位制:A=a[n]a[n-1]…a[0],B=b[n]b[n-1]…b[0]。
則組合數C(A,B)與C(a[n],b[n])C(a[n-1],b[n-1])

…*C(a[0],b[0]) mod p同餘
即:Lucas(n,m,p)=C(n%p,m%p)*Lucas(n/p,m/p,p)

#include<iostream>
//#include<algorithm>
using namespace std;
typedef long long ll;
int quick_power_mod(int a,int b,int m){//pow(a,b)%m
    int result = 1;
    int base = a;
    while(b>0){
         if(b & 1==1){
            result = (result*base) % m;
         }
         base = (base*base) %m;
         b>>=1;
    }
    return result;
}
//計算組合數取模
ll comp(ll a, ll b, int p) {//composite num C(a,b)%p
    if(a < b)   return 0;
    if(a == b)  return 1;
    if(b > a - b)   b = a - b;

    int ans = 1, ca = 1, cb = 1;
    for(ll i = 0; i < b; ++i) {
        ca = (ca * (a - i))%p;
        cb = (cb * (b - i))%p;
    }
    ans = (ca*quick_power_mod(cb, p - 2, p)) % p;
    return ans;
}
ll lucas(ll n, ll m, ll p) {
     ll ans = 1;

     while(n&&m&&ans) {
        ans = (ans*comp(n%p, m%p, p)) % p;//also can be recusive
        n /= p;
        m /= p;
     }
     return ans;
}
int main(){
    ll m,n;
    while(cin>>n>>m){
        cout<<lucas(n,m,10007)<<endl;
    }
    return 0;
}

C(n % mod, m % mod) % mod; 如果太大還是利用上面的逆元來處理。

半預處理
由於Lucas定理保證了階乘的數均小於p,所以可以講所有的階乘先預處理,優化C(n,m)
mod的要求:p<10^6,且為素數
有效範圍:1<=n,m<=10^9

//半預處理  
const ll MAXN = 100000;  
ll fac[MAXN+100];  
void init(int mod)  
{  
    fac[0] = 1;  
    for (int i=1; i<mod; i++){  
        fac[i] = fac[i-1] * i % mod;  
    }  
}  

//半預處理逆元求C(n,m)%mod  
ll C(ll n, ll m)  
{  
    if ( m>n ) return 0;  
    return fac[n] * (GetInverse(fac[m]*fac[n-m], mod)) % mod;  
}