Lucas定理與大組合數的取模的求法總結
分享一下我老師大神的人工智慧教程！零基礎，通俗易懂！http://blog.csdn.net/jiangjunshow
也歡迎大家轉載本篇文章。分享知識，造福人民，實現我們中華民族偉大復興！
Lucas定理與大組合數的取模的求法總結
分類： ACMer 數學 2012-03-11 09:38 1219人閱讀首先給出這個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以內解。
然後再大的資料，我就不會了！
貼一個例題的程式碼，方便以後看。
將不大於m顆種子存放在n顆樹中，問有多少種存法。
首先是不大於m顆種子，我沒可以認為少於m的那些種子存放在了第n+1顆樹上，這樣的話，問題就轉化成了將m顆種子存放在n+1顆樹上的方案數。ok這個是組合數學裡面的公式，亦即插板法，也就是X1+X2+X3+……+Xn+1 = m;ok，答案是C(n+m,m);
然後就是上面說的Lucas定理解決大組合數問題了
[cpp] view plain copy print ?
- #include <iostream>
- #include <string.h>
- #include <stdio.h>
- using namespace std;
- #define N 100010
- long long mod_pow(int a,int n,int p)
- {
- long long ret=1;
- long long A=a;
- while(n)
- {
- if (n & 1)
- ret=(ret*A)%p;
- A=(A*A)%p;
- n>>=1;
- }
- return ret;
- }
- long long factorial[N];
- void init(long long p)
- {
- factorial[0] = 1;
- for(int i = 1;i <= p;i++)
- factorial[i] = factorial[i-1]*i%p;
- //for(int i = 0;i < p;i++)
- //ni[i] = mod_pow(factorial[i],p-2,p);
- }
- long long Lucas(long long a,long long k,long long p) //求C(n,m)%p p最大為10^5。a,b可以很大！
- {
- long long re = 1;
- while(a && k)
- {
- long long aa = a%p;long long bb = k%p;
- if(aa < bb) return 0; //這個是最後的改動！
- re = re*factorial[aa]*mod_pow(factorial[bb]*factorial[aa-bb]%p,p-2,p)%p;//這兒的求逆不可先處理
- a /= p;
- k /= p;
- }
- return re;
- }
- int main()
- {
- int t;
- cin >> t;
- while(t--)
- {
- long long n,m,p;
- cin >> n >> m >> p;
- init(p);
- cout << Lucas(n+m,m,p) << "\n";
- }
- return 0;
- }
-
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定理推廣（組合數取模）
Lucas定理推廣（組合數取模）
2011年4月21日 3 條評論比賽的時候，意外想到了個方法來處理組合數取模。一推，發現就是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)的時間完成，我們就得到一個
時間複雜度的演算法（其中p是m的素因子，t是該素因子的個數），然後我們需要用lg級別的時間，完成每一個的合併，空間取決於素數表的大小（當然，可以沒有，沒有的話，就是lg級別的空間了）
可能的退化：如果p非常小，t非常大，即模數是小素數的高次。這個時候問題被退化到差不多O(k)的複雜度，k很大的情況下會很慢。但是經過測試，平均的情況非常快。
下面給出程式碼
CombinationNumber.cpp
包含一個比擴充套件lucas更穩定的版本（很長，包括了需要用到的所有模組，同餘方程組素數表什麼的都在，用的時候記得先生成素數表）：
PS:公式5： Wilson定理，當然這東西是什麼不重要，因為都會被約掉。
分類: ACM,Math Algorithm,數論,組合基於二進位制的遞推
2011年3月15日 沒有評論快半年沒更新了，囧。。
把去年的某類結果整理了下，總結如下：
我們都習慣了前後項的逐項遞推策略，這類一般的策略，可以稱之為線性方式。下面，我說一種當序列滿足某些性質是能狗構造的平方遞推模式(基於項的二進位制表示)
序列需要滿足的性質
***我們假設序列為a[1], a[2], a[3], a[4], a[5], ..., a[n]，以下所說的“簡單”，意思是計算時間代價低於線性
- a[n] -> a[2n]是簡單的（n是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]
- 假設p = n + m，n是不大於p的最大的2的冪，我們要做的只是求a[n]和a[m]。m大小的子問題和p大小的問題只在規模上存在不同。問題很快被遞迴到二進位制低位
- 現在我們倒過來，你懂的。按照二進位制展開逐位往高位遞推（而且，我們在遞推過程中逐位生成了a[n]，n是二的冪，真是一舉兩得）。
再給張圖片，方便大家理解，推導規模為1101的情況
時間效率分析
這個還是依賴於合併和地推的效率的，如果性質1和2的效率是O(1)的，那麼，最終效率就是O(lgn)，這個效率已經達到了資訊量上的線性，單從界的角度已經是最好的了。
例項
逐次平方求冪，移位乘法，等比前n項和等各種都可以用這種方法構造，效率就不用多說了，時間O(lgn)，空間O(1)
分類: ACM,Math組合數取模總結
2010年10月25日 8 條評論組合數取模，其實實際用途不大。。好吧，acm老喜歡出這樣的題目，因為資料可能很大。。
以下給出第一個約分版本的組合數取模，能夠處理當C(n,m) mod s當m不是很大的情況，有很多優化，先給程式碼吧^^
12345678910111213141516171819202122232425262728
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;}思路很直接——約分，基於這個公式
，直接吧分子全部儲存在一個數組裡，然後用下面的數不斷gcd，約去公因子，直到去約的數成為1，即下面所有的數都成為1。注意我加了一個優化——把分子和分母儘可能合併，以減少gcd執行的次數，這是一種貪心的策略。
基於以下的問題:給定若干容器，每個裡面存放了一個數，每個有一個容量的上限。規定一次合併操作——把兩個容器合併為一個，但是其中儲存的數為兩個數的乘積。求若干合併之後，容器總數最少。
策略：先排序，然後從最大的數開始，找最小的數和它乘，看有沒有溢位，沒有就合併這兩個數，並且“刪掉”最小的，用新的最小的數重複剛才的過程（和選定的那個大的數相乘），否則取次大的，不斷進行這個過程，直到到達序列末。不難證明它是最優的。基於這個優化，在m不是很大的時候（一般不超過100），這個演算法是接近線性(O(m))的，joj 2606就是這麼被我優化到0.00s的^^
第二個版本，使用素數分解約分，這個版本很快，估計應該是亞線性的，但它不能處理n非常大的情況，一般能夠處理n在最大1000000左右的情況，並且需要用到素數表。基於以下公式：
，先貼程式碼^^
123456789101112131415161718192021222324252627
///need prime table,and n < Nint 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的情況和模素數的乘方
的情況
12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667
///O(m) p is prime number, p*p < INT_MAXint 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_MAXint 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，則最後需要求的x值。
注意，如果b和p有大於1的公因子，會有多解，這個就是退化情況，所有在過程中需要不斷抽取b中p的因子，抽取之後再不斷取模，這樣就可以了。
當然，因為a中相應的要約去b中被抽取的因子(這些都被一些臨時變數儲存)，所以a也要抽取p的因子（需要注意的是，乘了再約和約了再乘是不一樣的，前者會導致失真）
當然，上面是通過同餘方程組，不用同餘方程組的話，可以得到一個m*lgn的演算法，不斷抽取模數m的因子^^，也給出程式碼吧，因為沒有同餘方程組，實現稍簡單一些，如果要求不是太高，也可以用
12345678910111213141516171819202122232425
/// s*s < INT_MAX , and s can be any positive numberint combination_Mod_s (int n,
首先給出這個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以內解。
然後再大的資料，我就不會了！
貼一個例題的程式碼，方便以後看。
將不大於m顆種子存放在n顆樹中，問有多少種存法。
首先是不大於m顆種子，我沒可以認為少於m的那些種子存放在了第n+1顆樹上，這樣的話，問題就轉化成了將m顆種子存放在n+1顆樹上的方案數。ok這個是組合數學裡面的公式，亦即插板法，也就是X1+X2+X3+……+Xn+1 = m;ok，答案是C(n+m,m);
然後就是上面說的Lucas定理解決大組合數問題了
[cpp] view plain copy print ?
- #include <iostream>
- #include <string.h>
- #include <stdio.h>
- using namespace std;
- #define N 100010
- long long mod_pow(int a,int n,int p)
- {
- long long ret=1;
- long long A=a;
- while(n)
- {
- if (n & 1)
- ret=(ret*A)%p;
- A=(A*A)%p;
- n>>=1;
- }
- return ret;
- }
- long long factorial[N];
- void init(long long p)
- {
- factorial[0] = 1;
- for(int i = 1;i <= p;i++)
- factorial[i] = factorial[i-1]*i%p;
- //for(int i = 0;i < p;i++)
- //ni[i] = mod_pow(factorial[i],p-2,p);
- }
- long long Lucas(long long a,long long k,long long p) //求C(n,m)%p p最大為10^5。a,b可以很大！
- {
- long long re = 1;
- while(a && k)