求组合数(取模)的两种方法

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

这里复述如下: 令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 − 1 , k − 1 ) + C ( n − 1 , k ) C ( n , 0 ) =   C ( n , n ) = 1 \left\{ \begin{aligned} C(n, k) = &\ C(n - 1, k - 1) + C(n - 1, k) \\ C(n, 0) = &\ C(n, n) = 1 \end{aligned} \right. {C(n,k)=C(n,0)= C(n1,k1)+C(n1,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;
    }
}

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

当然我们知道 C ( n , k ) = C ( n , n − k ) C(n, k) = C(n, n - k) C(n,k)=C(n,nk) ,所以上面的代码有很多空间和时间的浪费。可以将 tC[][] 二维数组转化为一维数组存储,同时,当 j &gt; i / 2 j &gt; i / 2 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);
    }
}

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

显然,由于空间的限制,pascal打表的方式并不适合求取一些比较大的组合数。例如,我们现在要求取的组合数的 n n n 的范围是 [ 1 , 1000000 ] [1, 1000 000] [1,1000000] , 那么我们应该怎么办呢? 那就轮到方法二大显身手了。
#法二:逆元求取组合数
由定理可知:如果用C(n, r)表示n-元素集的r-组合的个数,有 C ( n , r ) = n ! r ! ∗ ( n − r ) ! C(n, r) = \frac {n!} {r! * (n - r)!} C(n,r)=r!(nr)!n!

而我们的目标就是计算 C ( n , r ) % m o d C(n, r) \% mod 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 ≠ \neq ̸= (a % p / b % p) % p 。

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

首先给出数论中的解释:

对于正整数 a a a p p p,如果有 a x ≡ 1 ( m o d p ) ax \equiv 1 ( mod p) ax1(modp),那么把这个同余方程中 x x x 的最小正整数解叫做 a a a p p p 的逆元。

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

现在我们来解决模运算的除法问题。假设我们要计算的是 a b % p \frac a b \% p ba%p的值 m m m,即 a b % p = m      ( a 和 b 满 足 a % b = 0 ) \frac a b \% p = m \ \ \ \ (a和b满足 a \% b = 0) ba%p=m    (aba%b=0)

同时存在另一个数 x x x 满足 a x % p = m ax \% p = m ax%p=m

由模运算对乘法成立,两边同时乘以 b b b ,得到: a % p = ( m ( b % p ) ) % p a \% p = (m(b \% p)) \% p a%p=(m(b%p))%p

如果 a a a b b b 均小于模数 p p p 的话,上式可以改写为: a = b m % p a = bm \% p a=bm%p

等式两边再同时乘以 x x x, 得到: a x % p = m % p = x b m % p ax \% p = m \% p = xbm \% p ax%p=m%p=xbm%p

因此可以得到: b x % p = 1 bx \% p = 1 bx%p=1

哎,x是b的逆元呀(x 在模运算的乘法中等同于 1 b \frac 1 b b1, 这就是逆元的意义)

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

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

  1. 拓展欧几里得算法

  2. 费马小定理

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

对于第二种方法,因为在算法竞赛中模数p总是质数,所以可以利用费马小定理 b p − 1 % p = 1 b^{p-1} \% p = 1 bp1%p=1

可以直接得到 b b b 的逆元是 b p − 2 b^{p-2} bp2 , 使用 快速幂 求解即可。

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

  1. 求取1到n的阶乘对 mod 取模的结果存入数组 JC[] 中;
  2. 求取 C ( n , r ) C(n, r) C(n,r) 时, 先利用“拓展欧几里得算法”或者“费马小定理+快速幂”求 JC[r]的逆元存入临时变量 x 1 x_1 x1 ;
  3. 然后计算 J C [ n ] ∗ x 1 % m o d JC[n] * x_1 \% mod JC[n]x1%mod 存入临时变量 x 2 x_2 x2;( x 2 x_2 x2 即为 n ! r ! % m o d \frac {n!} {r!} \% mod r!n!%mod 的值)
  4. 求取JC[n - r] 的逆元存入临时变量 x 3 x_3 x3;
  5. 则可以得到 C ( n , r ) = x 2 ∗ x 3 % m o d C(n, r) = {x_2} * {x_3} \% mod C(n,r)=x2x3%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)
{
    if(a < b) return 0;
    return Jc[a] * niYuan(Jc[b], mod) % mod
        * niYuan(Jc[a - b], mod) % mod;
}

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

以上です。

评论 9
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值