Lucas定理与大组合数的取模的求法总结

 

Lucas定理与大组合数的取模的求法总结

分类: ACMer 数学 1219人阅读 评论(0) 收藏 举报
c

首先给出这个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以内解。

然后再大的数据,我就不会了!


贴一个例题的代码,方便以后看。

hdu 3037

将不大于m颗种子存放在n颗树中,问有多少种存法。

首先是不大于m颗种子,我没可以认为少于m的那些种子存放在了第n+1颗树上,这样的话,问题就转化成了将m颗种子存放在n+1颗树上的方案数。ok这个是组合数学里面的公式,亦即插板法,也就是X1+X2+X3+……+Xn+1 = m;ok,答案是C(n+m,m);

然后就是上面说的Lucas定理解决大组合数问题了

  1. #include <iostream>  
  2. #include <string.h>  
  3. #include <stdio.h>  
  4. using namespace std;  
  5.   
  6. #define N 100010  
  7.   
  8. long long mod_pow(int a,int n,int p)  
  9. {  
  10.     long long ret=1;  
  11.     long long A=a;  
  12.     while(n)  
  13.     {  
  14.         if (n & 1)  
  15.             ret=(ret*A)%p;  
  16.         A=(A*A)%p;  
  17.         n>>=1;  
  18.     }  
  19.     return ret;  
  20. }  
  21.   
  22. long long factorial[N];  
  23.   
  24. void init(long long p)  
  25. {  
  26.     factorial[0] = 1;  
  27.     for(int i = 1;i <= p;i++)  
  28.         factorial[i] = factorial[i-1]*i%p;  
  29.     //for(int i = 0;i < p;i++)  
  30.         //ni[i] = mod_pow(factorial[i],p-2,p);  
  31. }  
  32.   
  33. long long Lucas(long long a,long long k,long long p) //求C(n,m)%p p最大为10^5。a,b可以很大!  
  34. {  
  35.     long long re = 1;  
  36.     while(a && k)  
  37.     {  
  38.         long long aa = a%p;long long bb = k%p;  
  39.         if(aa < bb) return 0; //这个是最后的改动!  
  40.         re = re*factorial[aa]*mod_pow(factorial[bb]*factorial[aa-bb]%p,p-2,p)%p;//这儿的求逆不可先处理  
  41.         a /= p;  
  42.         k /= p;  
  43.     }  
  44.     return re;  
  45. }  
  46.   
  47. int main()  
  48. {  
  49.     int t;  
  50.     cin >> t;  
  51.     while(t--)  
  52.     {  
  53.         long long n,m,p;  
  54.         cin >> n >> m >> p;  
  55.         init(p);  
  56.         cout << Lucas(n+m,m,p) << "\n";  
  57.     }  
  58.     return 0;  
  59. }  
  60. Lucas定理的一个证明(找的)

    无意中看到这么一个定理,wiki上的证明我是看不懂了...

    http://en.wikipedia.org/wiki/Lucas%27_theorem 

    友情提示:上wiki,先翻墙- -||

    千辛万苦...终于找到一个我能看懂的证明~

    先贴个wiki上的公式:

     

    For non-negative integers m and n and a prime p, the following congruence relation holds:


    where


    and


    are the base p expansions of m and n respectively.

     

    首先我们注意到 n=(ak...a2,a1,a0)p  =  (ak...a2,a1)p * p + a0

                                                           =  [n/p]*p+a0


                                                        且m=[m/p]+b0

     

    只要我们更够证明 C(n,m)=C([n/p],[m/p]) * C(a0,b0)  (mod p)

    剩下的工作由归纳法即可完成

     

    我们知道对任意质数p:   (1+x)^p  == 1+(x^p)  (mod p) 

    注意!这里一定要是质数  ................(为什么)

     

    对 模p 而言

     上式左右两边的x^m的系数对模p而言一定同余(为什么),其中左边的x^m的系数是 C(n,m) 而由于a0和b0都小于p

    右边的x^m ( = x^(([m/p]*p)+b0)) 一定是由 x^([m/p]*p) 和 x^b0 相乘而得 (即发生于 i=[m/p] , j=b0 时) 因此我们就有了

     

    C(n,m)=C([n/p],[m/p]) * C(a0,b0)  (mod p) 

     

    这一结论!!!

     

    写的很乱~将就看吧

    Lucas定理推广(组合数取模)

    2011年4月21日Earthson3 条评论

    比赛的时候,意外想到了个方法来处理组合数取模。一推,发现就是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标签:Algorithm,数论,组合

    基于二进制的递推

    2011年3月15日Earthson没有评论

    快半年没更新了,囧。。

    把去年的某类结果整理了下,总结如下:

    我们都习惯了前后项的逐项递推策略,这类一般的策略,可以称之为线性方式。下面,我说一种当序列满足某些性质是能狗构造的平方递推模式(基于项的二进制表示)


    序列需要满足的性质

    ***我们假设序列为a[1], a[2], a[3], a[4], a[5], ..., a[n],以下所说的“简单”,意思是计算时间代价低于线性

    1. a[n] -> a[2n]是简单的(n是2的幂)
    2. a[n], a[m] -> a[n + m]是简单的(n是2的幂,且m小于n)

    当然,有一个更强的性质会使问题简化(它是上面条件的充分条件,但不是必要的)

    a[n], a[m] -> a[n + m]是简单的(m, n为任意正整数)


    算法构造

    首先考虑到因为性质1,所以,在已知a[1]推a[n]是非常简单的(其中n是2的幂)。
    其次,假设我们已经知道a[1]~a[n-1]的所有值,依赖性质2,我们可以轻易的把这些和a[n]合并,然后得到a[n + m]的值。于是,我们获得了a[1]~a[2n - 1]的所有值。
    注意到我们有a[1]所以,我们能对整个序列产生一个覆盖。

    接下去我们构造算法,来求特定的位置的a[p]

    1. 假设p = n + m,n是不大于p的最大的2的幂,我们要做的只是求a[n]和a[m]。m大小的子问题和p大小的问题只在规模上存在不同。问题很快被递归到二进制低位
    2. 现在我们倒过来,你懂的。按照二进制展开逐位往高位递推(而且,我们在递推过程中逐位生成了a[n],n是二的幂,真是一举两得)。

    再给张图片,方便大家理解,推导规模为1101的情况

    recurrence


    时间效率分析

    这个还是依赖于合并和地推的效率的,如果性质1和2的效率是O(1)的,那么,最终效率就是O(lgn),这个效率已经达到了信息量上的线性,单从界的角度已经是最好的了。


    实例

    逐次平方求幂,移位乘法,等比前n项和等各种都可以用这种方法构造,效率就不用多说了,时间O(lgn),空间O(1)

    分类:ACM,Math标签:

    组合数取模总结

    2010年10月25日Earthson8 条评论

    组合数取模,其实实际用途不大。。好吧,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标签:Algorithm,数论,组合

    素数筛法 Sieve Method

    2010年9月28日Earthson4 条评论

    前段时间学会了一个线性筛法,这个算法居然是O(n)的,确实令人震惊。而且更棒的是,它还为处理基于素因子分解的积性函数,提供了线性的筛法

    先贴一下代码吧^^

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    
    const int N=100000;
    int prime[N+1];
    int get_prime() ///get primes from 1 to N in liner time
    {
    	memset(prime,0,sizeof(int)*(N+1));
    	for(int i=2;i<=N;i++)
    	{
    		if(!prime[i]) prime[++prime[0]]=i;
    		for(int j=1;prime[j]<=N/i;j++)
    		{
    			prime[prime[j]*i]=1;
    			if(i%prime[j]==0) break;
    		}
    	}
    	return prime[0];
    }

    算法的思路是这样的:总是使用能整除合数的最小素因子,筛去合数。注意第12行,这句话保证i不被任何比prime[j]小的素数整除,对于i*prime[j](这个j包括比当前条件成立的j更小的j),它的最小素因子显然就是prime[j]。注意,筛去合数的,总是该合数的最小素因子和与之对应的另一个因子i,所以它只会被筛一次^^

    类似的可以得到很多积性以及“和性”的函数值,比如欧拉phi函数,sigma函数(n的所有因数和),还有n的所有因子的个数,等等等等。。。

    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
    
    int ph[N+1]={0};     	///phi(x) for x<=N
    int sig[N+1]={0};     	///sigma(x) for x<=N,sigma(x) is sum of x's all productors^^
    int productorcnt[N+1]; 	///count of productor numbers frome 1 ro N
    
    void phi_all() ///get phi(n) from 1 to N in liner time
    {
    	memset(ph,0,sizeof(int)*(N+1));
    	prime[0]=0;
    	for(int i=2;i<=N;i++)
    	{
    		if(!ph[i]) prime[++prime[0]]=i,ph[i]=i-1;
    		for(int j=1;prime[j]<=N/i;j++)
    		{
    			if(i%prime[j]) ph[prime[j]*i]=ph[prime[j]]*ph[i];
    			else
    			{
    			    ph[prime[j]*i]=ph[i]*prime[j];
    			    break;
    			}
    		}
    	}
    }
    
    
    int min_pk[N+1];
    void sigma_all()
    {
    	memset(sig,0,sizeof(int)*(N+1));
    	sig[1]=1,min_pk[1]=1,prime[0]=0;
    	for(int i=2;i<=N;i++)
    	{
    		if(!sig[i]) prime[++prime[0]]=i,sig[i]=i+1,min_pk[i]=i+1;
    		for(int j=1;prime[j]<=N/i;j++)
    		{
    			if(i%prime[j])
    			{
    			    sig[i*prime[j]]=sig[i]*sig[prime[j]];
    			    min_pk[i*prime[j]]=min_pk[prime[j]];
    			}
    			else
    			{
    				min_pk[i*prime[j]]=min_pk[i]*prime[j]+1;
    				sig[i*prime[j]]=(sig[i]/min_pk[i])*min_pk[i*prime[j]];
    				break;
    			}
    		}
    	}
    }

    注意:得到的是从1~N之间所有数的对应函数值,而且这上面的每一个函数,都能顺便产生一个对应素数表。用一个50000大小的表,你就能处理50000的平方大小的数的对应函数值了^^



没有更多推荐了,返回首页