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) ，特别的，Lucas（x，0，p）=1。

    其实说白了，Lucas定理就是求组合数C（n，m）mod p（p是素数）的值，即(（ n! /(n-m)!）/m! )mod p，而我们

又知道(a/b)mod p=a*b^(p-2)mod p（这里用到了一点逆元和费马小定理的知识），这样我们就可以在计算阶乘的过程中对p取模，不会造成溢出。（需要注意的是Lucas定理处理的p的范围大致为10^5数量级）

看下面一道题：HDU3037

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=3037

题目大意：m个种子要放在n颗树上，问有多少种方法，结果对素数p取模。

分析：m个种子可以分成两部分：放在树上的和不放在树上的，我们可以假想出第n+1颗树，认为那些没放在树上的种子放在这颗假想树上，这样，问题就变为了m个种子放在n颗树上有多少种方案数了。等价于方程X1+X2+...+Xn+1=m有多少组非负整数解，即(X1+1)+(X2+1)+...+(Xn+1+1)=m+n+1有多少组正整数解。挡板原理：（m+n+1）个1，分成n+1部分的方案数==>C（n+m,n）。到这儿就很明显了，果断Lucas。

实现代码如下：

#include <iostream>#include <cstdio>#include <cstring>using namespace std;typedef long long llg;const int N =150000;llg n, m, p, fac[N];void init(){    int i;    fac[0] =1;    for(i =1; i <= p; i++)        fac[i] = fac[i-1]*i % p;}llg quick_mod(llg a, llg b){    llg tmp = a % p, ans =1;    while(b)    {        if(b &1)  ans = ans * tmp % p;        tmp = tmp*tmp % p;        b >>=1;    }    return  ans;}llg C(llg n, llg m){    if(m > n)  return 0;    return  fac[n]*quick_mod(fac[m]*fac[n-m], p-2) % p;}llg Lucas(llg n, llg m){    if(m ==0)  return 1;    else return  (C(n%p, m%p)*Lucas(n/p, m/p))%p;}int main(){    int t;    scanf("%d", &t);    while(t--)    {        scanf("%I64d%I64d%I64d", &n, &m, &p);        init();        printf("%I64d\n", Lucas(n+m, m));    }    return 0;}

对于p比较大的情况，不能对阶乘预先处理，需要单独处理。

题目链接：http://acm.fzu.edu.cn/problem.php?pid=2020

题目大意：求C（n，m）mod p。n,p是10^9数量级的。

实现代码如下：

#include <iostream>#include <string.h>#include <stdio.h>using namespace std;typedef long long LL;LL n,m,p;LL quick_mod(LL a, LL b){    LL ans = 1;    a %= p;    while(b)    {        if(b & 1)        {            ans = ans * a % p;            b--;        }        b >>= 1;        a = a * a % p;    }    return ans;}LL C(LL n, LL m){    if(m > n) return 0;    LL ans = 1;    for(int i=1; i<=m; i++)    {        LL a = (n + i - m) % p;        LL b = i % p;        ans = ans * (a * quick_mod(b, p-2) % p) % p;    }    return ans;}LL Lucas(LL n, LL m){    if(m == 0) return 1;    return C(n % p, m % p) * Lucas(n / p, m / p) % p;}int main(){    int T;    scanf("%d", &T);    while(T--)    {        scanf("%I64d%I64d%I64d", &n, &m, &p);        printf("%I64d\n", Lucas(n,m));    }    return 0;}