存档

文章标签 ‘组合’

Lucas定理推广（组合数取模）

2011年4月21日 7 条评论

比赛的时候，意外想到了个方法来处理组合数取模。一推，发现就是Lucas定理，并且，在模数是素数乘方的情况作了些推广。使得lucas能够处理模任意数的组合数取模。

Lucas定理

首先给出Lucas定理的递归描述

lucas定理

简要证明思路：

考虑连续的p个整数（p是素数），把能整除p的那个数除去，然后取模，再乘起来模p其实是-1(*)。注意，因为和p互素的原因，这东西是可约的。

再考虑这个组合公式

上下都有k项，你能取的长度为p的连续整数有k/p段(从后往前取)，每一段去掉那个能被p整除的数，然后，上下因为有相同的段数，所以就被约掉了（明白为什么在这种情况下是可约的么？）。

为什么能被p整除的数要被单独的拉出来呢？因为，和模数有公因子的数加入取模运算会导致不可逆。所以这些能被p整除的数需要被提取出来单独处理。

于是，有了后面部分的C(n/p, k/p)，其实就是每个被提取出来的数再提取了一个因子p（因为上下提取出来的数一样多，所以，这些p自然就被约掉了）。然后，你幸运的发现，他是个组合数，于是问题被递归。最后，我们还有个前缀没有被处理。但是，因为模数是p，所以有下面公式

把这写合并在一起，就有了Lucas定理的递归表述（公式1）

推广的Lucas

那么要怎么处理模数是素数p的乘方的情形呢？其实是一样的，只不过，有些东西就不再退化了，比如公式6。我们可以用同样的方式，先把因子p全部提取出来单独处理，于是得到如下的公式，来处理p的乘方。

其中

式2的前缀（F除以F）就是上面的公式6不能退化的情况。为什么？因为其中可能含有因子p。因子p加入会导致乘法取模的不可逆。

需要注意的是，前缀很可能不是整数，在写代码的时候这个地方很容易犯错！所以要先进入子问题（可以用一个栈搞定）。

一般的组合数取模

有了这些，我们能干什么呢？其实，Lucas及其扩展只是把原来的组合数取模问题分割成了更小规模的问题（有lg(k)/lg(p)个）。

考虑一个一般的问题，模数m任意。

策略：

我们把m先素数分解，然后用Lucas及其扩展分割。然后对每个小问题解决合并。合并可以用CRT进行（同余方程组，中国剩余定理）。于是，如果简单的组合数取模能在O(k)的时间完成，我们就得到一个O(\min (k, p^t \cdot log_p(t)))时间复杂度的算法（其中p是m的素因子，t是该素因子的个数），然后我们需要用lg级别的时间，完成每一个的合并，空间取决于素数表的大小（当然，可以没有，没有的话，就是lg级别的空间了）

可能的退化：如果p非常小，t非常大，即模数是小素数的高次。这个时候问题被退化到差不多O(k)的复杂度，k很大的情况下会很慢。但是经过测试，平均的情况非常快。

下面给出代码
CombinationNumber.cpp
包含一个比扩展lucas更稳定的版本

（很长，包括了需要用到的所有模块，同余方程组素数表什么的都在，用的时候记得先生成素数表）：

PS:公式5： Wilson定理，当然这东西是什么不重要，因为都会被约掉。

分类: ACM, Math 标签: , ,

组合数取模总结

2010年10月25日 8 条评论

组合数取模，其实实际用途不大。。好吧，acm老喜欢出这样的题目，因为数据可能很大。。

以下给出第一个约分版本的组合数取模，能够处理当C(n,m) mod s当m不是很大的情况，有很多优化，先给代码吧^^

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

int combination_Mod_m (int a, int b, int m)
{
    const int gate = 2000000000;  ///larger gate for more Optimization
    if (b > a) return 0;
    b = (a - b < b) ? a - b : b;
    int d[b], tmp = 1 ,l = 0, h, j, g;
    for (l = 0, j = a - b + 1; a - l >= j; l++)
    {
		d[l] = a - l;
		while(gate / d[l] >= j && a - l != j)
			d[l] *= j, j++;
    }
    for (j = 2, h = b; j <= h; h--)
    {
    	tmp *= h;
        while (gate / tmp >= j && h != j)
        	tmp *= j, j++;
        for (int i = 0; tmp != 1; i++)
        {
            g = gcd (tmp, d[i]);
            d[i] /= g, tmp /= g;
        }
    }
    int result = 1;
    for (int i = 0; i < l; i++)
        d[i] %= m, result *= d[i], result %= m;
    return result;
}

思路很直接——约分，基于这个公式(^n _m) = \frac{n \cdot (n-1) \cdots (n-m+1)}{1 \cdot 2 \cdot (m-1)\cdots m}，直接吧分子全部保存在一个数组里，然后用下面的数不断gcd，约去公因子，直到去约的数成为1，即下面所有的数都成为1。注意我加了一个优化——把分子和分母尽可能合并，以减少gcd运行的次数，这是一种贪心的策略。

基于以下的问题:给定若干容器，每个里面存放了一个数，每个有一个容量的上限。规定一次合并操作——把两个容器合并为一个，但是其中保存的数为两个数的乘积。求若干合并之后，容器总数最少。

策略：先排序，然后从最大的数开始，找最小的数和它乘，看有没有溢出，没有就合并这两个数，并且“删掉”最小的，用新的最小的数重复刚才的过程（和选定的那个大的数相乘），否则取次大的，不断进行这个过程，直到到达序列末。不难证明它是最优的。基于这个优化，在m不是很大的时候（一般不超过100），这个算法是接近线性(O(m))的，joj 2606就是这么被我优化到0.00s的^^

第二个版本，使用素数分解约分，这个版本很快，估计应该是亚线性的，但它不能处理n非常大的情况，一般能够处理n在最大1000000左右的情况，并且需要用到素数表。基于以下公式：(^n _m) = \frac{n!}{m!\cdot(n-m)!}，先贴代码^^

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

///need prime table,and n < N
int combination_Mod_h (int n, int m, int h)
{
	if(h == 1) return 0;
        int result = 1, cnt = 0, temp;
	for(int i = 1; i < prime[0] && prime[i] <= n; i++)
	{
		temp = n, cnt = 0;
		while(temp)
			temp /= prime[i], cnt += temp;
		temp = n - m;
		while(temp)
			temp /= prime[i], cnt -= temp;
		temp = m;
		while(temp)
			temp /= prime[i], cnt -= temp;
		temp = prime[i];
		while(cnt)
		{
			if(cnt & 1)
                                result *= temp, result %= h;
			temp *= temp, cnt >>= 1, temp %= h;
		}
		if(result == 0) return 0;
	}
	return result;
}

这个基本没啥好说的，都比较直接，有一点要提一下：m!中素因子p的数目，求法:m不断除p，除的结果不断加到一个变量保存，当m为0的时候，那个保存的结果即所求（第9,12,15行都是这个过程），如果你大脑能作一个数域的映射的话，这个是显然的，当然你也可以技巧性的通过公式变形得到，具体就不解释了^^

接下去介绍一个线性O(m)的方法，简单的说，就是把模数分解，然后同余方程组合并，这里给出模素数p的情况和模素数的乘方p^k的情况

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
54
55
56
57
58
59
60
61
62
63
64
65
66
67

///O(m) p is prime number, p*p < INT_MAX
int combination_Mod_p (int n, int m, int p)
{
    if (m > n) return 0;
    m = (n - m < m) ? n - m : m;
    int a = 1, b = 1, x, y, pcnt = 0;
    for (int i = 1; i <= m; i++)
    {
        a *= n - i + 1, b *= i;
        while (a % p == 0) a /= p, pcnt++;
        while (b % p == 0) b /= p, pcnt--;
        b %= p, a %= p;
    }
	if (pcnt) return 0;
    extended_gcd(b, p, x, y);
    if (x < 0) x += p;
    x *= a, x %= p;
    return x;
}
///O(m)
///p is prime and p^k < INT_MAX
int combination_Mod_pk (int n, int m, int p, int k)
{
    if (m > n) return 0;
    m = (n - m < m) ? n - m : m;
    int a = 1, b = 1, x, y, pa = 1, pb = 1, pcnt = 0;
    int q = power(p, k);
    for (int i = 1; i <= m; i++)
    {
        a *= n - i + 1, b *= i;
        while(a % p == 0) pa *= p, a /= p;
        while(b % p == 0) pb *= p, b /= p;
        while (pa % q == 0) pa /= q, pcnt++;
        while (pb % q == 0) pb /= q, pcnt--;
        b %= q, a %= q;
    }
    if(pa < pb) pcnt--, pa *= q;
    pa /= pb, a *= pa, a %= q;
    if (pcnt) return 0;
    extended_gcd(b, q, x, y);
    if (x < 0) x += q;
    x *= a, x %= q;
    return x;
}
/**
  * Extended Euclid's Algorithm
  * ax+by=gcd(a,b);
  **/
int extended_gcd (int a, int b, int &x, int &y)
{
    if (!b)
        return x = 1, y = 0, a;
    int ret = extended_gcd (b, a % b, x, y), tmp = x;
    x = y, y = tmp - (a / b) * y;
    return ret;
}
int power (int x, int k)
{
    int result = 1;
    while (k)
    {
        if (k & 1)
            result *= x;
        x *= x, k >>= 1;
    }
    return result;
}

这个原理不难，就是分子分母分别对p取模，然后求逆。。当然，需要注意可能的退化情况。
假设分母为b，分子是a，则最后需要求a \equiv b \cdot x\, (mod \, p)的x值。
注意，如果b和p有大于1的公因子，会有多解，这个就是退化情况，所有在过程中需要不断抽取b中p的因子，抽取之后再不断取模，这样就可以了。
当然，因为a中相应的要约去b中被抽取的因子(这些都被一些临时变量保存)，所以a也要抽取p的因子（需要注意的是，乘了再约和约了再乘是不一样的，前者会导致失真）

当然，上面是通过同余方程组，不用同余方程组的话，可以得到一个m*lgn的算法，不断抽取模数m的因子^^，也给出代码吧，因为没有同余方程组，实现稍简单一些，如果要求不是太高，也可以用

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

/// s*s < INT_MAX , and s can be any positive number
int combination_Mod_s (int n, int m, int s)
{
    if (m > n) return 0;
    m = (n - m < m) ? n - m : m;
    int a = 1, b = 1, x, y, g, pcnt = 0;
    int pa = 1, pb = 1; ///To be sure that pa and pb is not overflow
    for (int i = 1; i <= m; i++)
    {
        a *= n - i + 1, b *= i;
        while( (g = gcd(a, s)) > 1) pa *= g, a /= g;
        while( (g = gcd(b, s)) > 1) pb *= g, b /= g;
        g = gcd(pa, pb), pa /= g, pb /= g;
        while (pa % s == 0) pa /= s, pcnt++;
        while (pb % s == 0) pb /= s, pcnt--;
        b %= s, a %= s;
    }
    while(pcnt) pcnt--, pa *= s;
    pa /= pb, pa %= s, a *= pa, a %= s;
    if (pcnt) return 0;
    extended_gcd(b, s, x, y);
    if (x < 0) x += s;
    x *= a, x %= s;
    return x;
}

最后，是一个我还没完全参透的定理：Lucas定理，具体可以参看wiki，目前还不太清楚能不能处理p的乘方（要是能就好了^^），慢慢研究吧

1
2
3
4
5
6
7
8
9
10
11

/**
  * Lucas' theorem---combination number mod prime p
  * Be sure that p*p<INT_MAX
  **/
int lucas (int m, int n, int p)
{
    int result = 1;
    while (m && n && result)
        result *= combination_Mod_p (m % p, n % p, p), result %= p, m /= p, n /= p;
    return result;
}
分类: ACM, Math 标签: , ,

组合数模p的计算

2010年4月6日 没有评论

这个计算的目的是什么？当组合数较大的时候，显然不能很好保存，因此把它转化为p的模（p是任意的）。

方案很简单，仍然利用公式的除法优先（我的另一篇文章组合数计算），不同的是，这次需要保证除法的绝对优先，这将导致不那么高的效率

首先明确一点a*b/c mod t和(a mod t)*b/c mod t可能是不同的，因为如果a能被c整除，a mod t却不一定能被c整除，(a mod t)*b/c这个式子甚至的结果甚至不是整数（如果那个除号不是整除的话。）比如：4*5/2 mod 3=1，而 4 mod 3=1 ，1*5/2不是一个整数。如果你看了我的那篇文章就知道，在那个算法中我对这一点的依赖

因此，我选择用数组保存了公式上边的因子。现在，我的目标是，把下面的因子全部约掉，我使用gcd（最大公约数）找上面元素和下面元素的因子，并且约掉该因子，扫描……显然最后下面的因子全部归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
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

int gcd(int a,int b)
{
    int swap;
    while(b!=0)
    {
        swap=b;
        b=a%b;
        a=swap;
    }
    return a;
}
 
int combinationNumMod_p(int n,int k,int p)
{
    int solution=1;
    int *dividend;
    k=(n-k>k)?k:n-k;
    dividend=new int[k];
    for(int i=0;i<k;i++)
    {
        dividend[i]=n-i;
    }
    for(int i=2;i<=k;i++)   //i就是生成的下因子
    {
        int s,g;
        s=i;                 //把i的值转给s，因i的值不能随便改（i负责生成下一个s）
        for(int j=0;j<k;j++)    //对上因子扫描
        {
            if(s==1)               //如果当前下因子已经是1了，就跳出循环
            {
                break;
            }
            g=gcd(s,dividend[j]); //这个算法的还依赖于gcd，g是保存的最大公约数
            if(g!=1)       //公约数不等于1就把公约数约掉
            {
                s/=g;                
                dividend[j]/=g;
            }
        }
    }
    for(int i=0;i<k;i++)        //把处理好的上因子乘起来
    {
        solution*=dividend[i]%p;      //公式a*b mod p=(a mod p)*(b mod p) mod p
        solution%=p;
    }
    return solution;
}
分类: Cpp, Math 标签:

可重复的组合

2010年4月3日 没有评论

这是一个困扰我很久的问题，刚上高中就开始找它的计数方案（证明），直到高三才搞定。我的把这种计数的情况称为潜在空位和潜在禁止位，它能解决一类组合计数问题。

当然，先看一下，可重复的组合。

有n种元素的一个集合，每种的数目假设无限（至少比k大），现在我们要选取其中的k个元素，令该组合所有情况的集合为T，求组合数。在书上我们很容易查到，它等于C(n+k-1,k)

先考虑一种一般的策略（在非可重复组合中）：有n个座位，我要选择其中的k个（我把这个集合称为C(n,k)），为了方便，同时我把座位标号（1，2，3，……），我可以这样做，首先选择一个座位，然后在该座的后边（标号更大）选择第二个，……依次选择，直到最后一个。很容易证明它没有重复，它对C(n,k)能建立一一映射

对于可重复组合，我只是加入了一个潜在位而已，比如，当我选择1号位，那么在1号位的后面将产生一个新的1号位，可供后面选择，选择k位的情况将会产生一个n+k个位子的现象，但是请注意，选择的第k个位子后面产生的潜在位，不能供任何其它位选择（或者说，那个位子在任何情况下都不会被选择到），所以实际上该位不是一个供选择的位。这样产生了k-1个“虚位”，很容易建立它与C(n+k-1,k)的一一映射。而且，它对T也是一一映射的（在下图中，最后一句话的意思是，在n+k-1位的等价图中，不存在两种不同情况，他们解释到T的时候，是一样的）

现在来看另一种情况，即我所谓的潜在禁止位：当我选择位子的时候，有这样的要求，相邻位子的必须有一个最小间隔，甚至不同位子间最小间隔可以不同。这其实并不困难，我们可以这样理解：每个选择了的位子后面有一个给定长度s的禁止位子，当我选择某个位子的时候，后面s个位子就等于消失了。（或者你可以理解为我选择的位子长度不是1，而是s+1……）这与重复组合的情况是类似的，同样注意最后一个被选择的位子。容易得到，对于给定的最小间隔s，有公式C(n-s(k-1),k),其中的s(k-1)是对禁止位的计数（很奇怪吧，这些禁止位的位置是不确定的，但他们确实存在，而且就是有那么多位，对于不同间隔的情况，把它的禁止位数替换上去就可以了）

分类: Math 标签:

组合数的计算

2010年4月2日 2 条评论

其实我们都知道组合有一个公式：

可以用于递归计算，但是，事实告诉我，这不是一个好算法，效率不行，因此需要一个更高效的算法。

显然注意到公式

它的效率显然太高了，悲剧的是如果用上面的乘积除以下面的乘积，很快就溢出了！这样做显然不行。

怎样才能避免溢出呢？对策：让除法优先。

那么怎么保证作用除法的时候一定能整除呢？注意到以下事实，连续的2个整数乘积，必然能被2!整除，连续的3个整数乘积，必然能被3!整除……连续k个整数乘积，必然能被k!整除。其实假设设连续的那k个整数最大的为n，那么整除结果必然是C(n,k)，它显然是一个整数.而且，只要最终结果不溢出，它就不会溢出！

于是给出 代码：

int combinationNumber(int n,int m)
{
    int solution=1,g;
    m=(n-m<m)?n-m:m; //一个优化，基于C(n,m)=C(n,n-m)
    for(int i=0;i<m;i++)
    {
        if(solution%(i+1)==0)
        {
            solution/=i+1;
            solution*=n-i;
        }
        else if((n-i)%(i+1)==0)
        {
            solution*=(n-i)/(i+1);
        }
        else   //之前没考虑到的情况，修复了这个bug
        {
            g=gcd(solution,i+1);
            solution/=g;
            solution*=(n-i)/((i+1)/g);
        }
    }
    return solution;
}
分类: Cpp 标签: