Lucas定理与大组合数的取模

首先给出这个Lucas定理：

A、B是非负整数，p是质数。AB写成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])  modp同余

即：Lucas(n,m,p)=c(n%p,m%p)*Lucas(n/p,m/p,p) 

这个定理的证明不是很简单，我一直想找个很好的张明，但是，没找到，昨天看到了一个解题报告，基本上可以说明白这个Lucas定理是怎么回事了，具体的是说：

以求解n! % p为例，把n分段，每p个一段，每一段求的结果是一样的。但是需要单独处理每一段的末尾p, 2p, ...，把p提取出来，会发现剩下的数正好又是(n / p)!，相当于划归成了一个子问题，这样递归求解即可。

这个是单独处理n!的情况，当然C(n,m)就是n!/(m!*(n-m)!)，每一个阶乘都用上面的方法处理的话，就是Lucas定理了，注意这儿的p是素数是有必要的。

Lucas最大的数据处理能力是p在10^5左右，不能再大了，hdu 3037就是10^5级别的！

对于大组合数取模，n,m不大于10^5的话，用逆元的方法，可以解决。对于n,m大于10^5的话，那么要求p<10^5，这样就是Lucas定理了，将n,m转化到10^5以内解。

概述

首先我们要知道什么是组合数。具体可以参考我之前的博客 “排列与组合”笔记 中，集合的组合的部分。

这里复述如下： 令r为非负整数。我们把n个元素的集合S的r-组合理解为从S的n个元素中对r个元素的无序选择。换句话说，S的一个r-组合是S的一个子集，该子集由S的n个元素中的r个组成，即S的一个r-元素子集。

由此，求解组合数即变成了求式子C(n, r) 的值。

法一：Pascal公式打表

由Pascal公式（参考 组合数学笔记之二——二项式系数），我们知道 

{C(n,k)=C(n,0)= C(n−1,k−1)+C(n−1,k) C(n,n)=1

取二维数组 tC[][] ，初始化 tC[0][0] = 1; 打表即可。代码最简单，如下：

const int maxn(1005), mod(100003);
int tC[maxn][maxn]; //tC 表示 table of C

inline int C(int n, int k)
{
    if(k > n) return 0;
    return tC[n][k];
}

void calcC(int n)
{
    for(int i = 0; i < n; i++)
    {
        tC[i][0] = 1;
        for(int j = 1; j < i; j++)
            tC[i][j] = (C(i - 1, j - 1) + C(i - 1, j)) % mod;
        tC[i][i] = 1;
    }
}

 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

计算  C(n,k)  返回内联函数 C(n,k)  的值即可。

当然我们知道  C(n,k)=C(n,n−k)  ，所以上面的代码有很多空间和时间的浪费。可以将 tC[][] 二维数组转化为一维数组存储，同时，当  j>i/2  时终止第二层循环，新代码如下：

const int maxn(10005), mod(100003);
int tC[maxn * maxn]; //tC 表示 table of C

inline int loc(int n, int k) // C(n, k)返回在一维数组中的位置
{
    int locate = (1 + (n >> 1)) * (n >> 1); // (n >> 1) 等价于 (n / 2)
    locate += k;
    locate += (n & 1) ? (n + 1) >> 1 : 0; // (n & 1) 判断n是否为奇数
    return locate;
}

inline int C(int n, int k)
{
    if(k > n) return 0;
    k = min(n - k, k);
    return tC[loc(n, k)];
}

void calcC(int n)
{
    for(int i = 0; i < n; i++)
    {
        tC[loc(i, 0)] = 1;
        for(int j = 1, e = i >> 1; j <= e; j++)
            tC[loc(i, j)] = C(i - 1, j) + C(i - 1, j - 1);
    }
}

 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

同样，要得到  C(n,k)  只需要返回内联函数 C(n,k)  的值即可。

显然，由于空间的限制，pascal打表的方式并不适合求取一些比较大的组合数。例如，我们现在要求取的组合数的  n  的范围是  [1,1000000]  , 那么我们应该怎么办呢？ 那就轮到方法二大显身手了。

法二：逆元求取组合数

由定理可知：如果用C(n, r)表示n-元素集的r-组合的个数，有

C(n,r)=n!r!∗(n−r)!

而我们的目标就是计算  C(n,r)%mod  的值。

由数论的知识我们知道，模运算的加法，减法，乘法和四则运算类似，即： 
模运算与基本四则运算有些相似，但是除法例外。其规则如下：

• (a + b) % p = (a % p + b % p) % p
• (a - b) % p = (a % p - b % p) % p
• (a * b) % p = (a % p * b % p) % p

但对于除法却不成立，即(a / b) % p  ≠  (a % p / b % p) % p 。

显然数学家们是不能忍受这种局面的，他们扔出了“逆元”来解决这个问题。那么什么是逆元？ 逆元和模运算中的除法又有说明关系呢？

首先给出数论中的解释：

对于正整数  a  和  p ，如果有  ax≡1(modp) ，那么把这个同余方程中  x  的最小正整数解叫做  a  模  p  的逆元。

什么意思呢？ 就是指，如果  ax%p=1  ， 那么x的最小正整数解就是  a  的逆元。

现在我们来解决模运算的除法问题。假设

ab

同时存在另一个数  x  满足

ax%p=m

由模运算对乘法成立，两边同时乘以  b  ，得到：

a%p=(m(b%p))%p

如果  a  和  b  均小于模数  p  的话，上式可以改写为：

a=bm%p

等式两边再同时乘以  x , 得到：

ax%p=m%p=xbm%p

因此可以得到：

bx%p=1

哎，x是b的逆元呀（x 在模运算的乘法中等同于  1b , 这就是逆元的意义）

由以上过程我们看到，求取  (ab%p)  等同于 求取  (a∗(b的逆元)%p)  。 因此，求模运算的除法问题就转化为就一个数的逆元问题了。

而求取一个数的逆元，有两种方法

1. 拓展欧几里得算法

2. 费马小定理

对于利用拓展欧几里得算法求逆元，很显然，如果 bx%p=1 ，那么  bx+py=1 , 直接利用 exgcd(b, p, x, y)（代码实现在后面给出），则  (x%p+p)%p  即为  b  的逆元。

对于第二种方法，因为在算法竞赛中模数p总是质数，所以可以利用费马小定理 ：

bp−1%p=1

可以直接得到  b  的逆元是  bp−2  , 使用 快速幂 求解即可。

明白了以上几个关键点，那么求取组合数  C(n,r)  的算法就呼之欲出了：

1. 求取1到n的阶乘对 mod 取模的结果存入数组 JC[] 中；
2. 求取  C(n,r)  时， 先利用“拓展欧几里得算法”或者“费马小定理+快速幂”求 JC[r]的逆元存入临时变量  x1  ;
3. 然后计算 JC[n]∗x1%mod  存入临时变量  x2 ;（ x2  即为 n!r!%mod  的值）
4. 求取JC[n - r] 的逆元存入临时变量  x3 ;
5. 则可以得到  C(n,r)=x2∗x3%mod

下面是方法二的代码片段：

typedef long long LL;
const LL maxn(1000005), mod(1e9 + 7);
LL Jc[maxn];

void calJc()    //求maxn以内的数的阶乘
{
    Jc[0] = Jc[1] = 1;
    for(LL i = 2; i < maxn; i++)
        Jc[i] = Jc[i - 1] * i % mod;
}
/*
//拓展欧几里得算法求逆元
void exgcd(LL a, LL b, LL &x, LL &y)    //拓展欧几里得算法
{
    if(!b) x = 1, y = 0;
    else
    {
        exgcd(b, a % b, y, x);
        y -= x * (a / b);
    }
}

LL niYuan(LL a, LL b)   //求a对b取模的逆元
{
    LL x, y;
    exgcd(a, b, x, y);
    return (x + b) % b;
}
*/

//费马小定理求逆元
LL pow(LL a, LL n, LL p)    //快速幂 a^n % p
{
    LL ans = 1;
    while(n)
    {
        if(n & 1) ans = ans * a % p;
        a = a * a % p;
        n >>= 1;
    }
    return ans;
}

LL niYuan(LL a, LL b)   //费马小定理求逆元
{
    return pow(a, b - 2, b);
}

LL C(LL a, LL b)    //计算C(a, b)
{
    return Jc[a] * niYuan(Jc[b], mod) % mod
        * niYuan(Jc[a - b], mod) % mod;
}
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

以上即为逆元求取组合数的方法，无论使用拓展欧几里得还是费马小定理，一开始求取Jc数组是的复杂度是  O(n) ，拓展欧几里得算法和费马小定理的复杂度均为  O(lg(mod))  , 如果要求取m次组合数，则总的时间复杂度为  O(mnlg(mod))