[关闭]
@ArrowLLL 2017-04-06T10:55:58.000000Z 字数 3035 阅读 1612

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

数学 组合数学 算法


主页地址 :月光森林


概述

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

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

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

法一:Pascal公式打表

由Pascal公式(参考组合数学笔记之二——二项式系数),我们知道


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

  1. const int maxn(1005), mod(100003);
  2. int tC[maxn][maxn]; //tC 表示 table of C
  3. inline int C(int n, int k)
  4. {
  5. if(k > n) return 0;
  6. return tC[n][k];
  7. }
  8. void calcC(int n)
  9. {
  10. for(int i = 0; i < n; i++)
  11. {
  12. tC[i][0] = 1;
  13. for(int j = 1; j < i; j++)
  14. tC[i][j] = (C(i - 1, j - 1) + C(i - 1, j)) % mod;
  15. tC[i][i] = 1;
  16. }
  17. }

计算 返回内联函数 的值即可。

当然我们知道 ,所以上面的代码有很多空间和时间的浪费。可以将 tC[][] 二维数组转化为一维数组存储,同时,当 时终止第二层循环,新代码如下:

  1. const int maxn(10005), mod(100003);
  2. int tC[maxn * maxn]; //tC 表示 table of C
  3. inline int loc(int n, int k) // C(n, k)返回在一维数组中的位置
  4. {
  5. int locate = (1 + (n >> 1)) * (n >> 1); // (n >> 1) 等价于 (n / 2)
  6. locate += k;
  7. locate += (n & 1) ? (n + 1) >> 1 : 0; // (n & 1) 判断n是否为奇数
  8. return locate;
  9. }
  10. inline int C(int n, int k)
  11. {
  12. if(k > n) return 0;
  13. k = min(n - k, k);
  14. return tC[loc(n, k)];
  15. }
  16. void calcC(int n)
  17. {
  18. for(int i = 0; i < n; i++)
  19. {
  20. tC[loc(i, 0)] = 1;
  21. for(int j = 1, e = i >> 1; j <= e; j++)
  22. tC[loc(i, j)] = C(i - 1, j) + C(i - 1, j - 1);
  23. }
  24. }

同样,要得到 只需要返回内联函数 的值即可。

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

法二:逆元求取组合数

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

而我们的目标就是计算 的值。

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

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

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

首先给出数论中的解释:

对于正整数 ,如果有 ,那么把这个同余方程中 的最小正整数解叫做 的逆元。

什么意思呢? 就是指,如果 , 那么x的最小正整数解就是 的逆元。

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

同时存在另一个数 满足

由模运算对乘法成立,两边同时乘以 ,得到:

如果 均小于模数 的话,上式可以改写为:

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

因此可以得到:

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

由以上过程我们看到,求取 等同于 求取 。 因此,求模运算的除法问题就转化为就一个数的逆元问题了。

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

  1. 拓展欧几里得算法

  2. 费马小定理

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

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

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

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

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

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

  1. typedef long long LL;
  2. const LL maxn(1000005), mod(1e9 + 7);
  3. LL Jc[maxn];
  4. void calJc() //求maxn以内的数的阶乘
  5. {
  6. Jc[0] = Jc[1] = 1;
  7. for(LL i = 2; i < maxn; i++)
  8. Jc[i] = Jc[i - 1] * i % mod;
  9. }
  10. /*
  11. //拓展欧几里得算法求逆元
  12. void exgcd(LL a, LL b, LL &x, LL &y) //拓展欧几里得算法
  13. {
  14. if(!b) x = 1, y = 0;
  15. else
  16. {
  17. exgcd(b, a % b, y, x);
  18. y -= x * (a / b);
  19. }
  20. }
  21. LL niYuan(LL a, LL b) //求a对b取模的逆元
  22. {
  23. LL x, y;
  24. exgcd(a, b, x, y);
  25. return (x + b) % b;
  26. }
  27. */
  28. //费马小定理求逆元
  29. LL pow(LL a, LL n, LL p) //快速幂 a^n % p
  30. {
  31. LL ans = 1;
  32. while(n)
  33. {
  34. if(n & 1) ans = ans * a % p;
  35. a = a * a % p;
  36. n >>= 1;
  37. }
  38. return ans;
  39. }
  40. LL niYuan(LL a, LL b) //费马小定理求逆元
  41. {
  42. return pow(a, b - 2, b);
  43. }
  44. LL C(LL a, LL b) //计算C(a, b)
  45. {
  46. return Jc[a] * niYuan(Jc[b], mod) % mod
  47. * niYuan(Jc[a - b], mod) % mod;
  48. }

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

以上です~

添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注