ACM數論之旅10---大組合數-盧卡斯定理（在下盧卡斯，你是個人Master嗎？(。-`ω´-) ）

 

記得前幾章的組合數吧

咱們學了O(n^2)的作法，加上逆元，咱們又會了O(n)的作法

如今來了新問題，若是n和m很大呢，

好比求C(n, m) % p  ， n<=1e18,m<=1e18,p<=1e5

看到沒有，n和m這麼大，可是p卻很小，咱們要利用這個p

 

（數論就是這麼無聊的東西，我要是讓n=1e100，m=1e100，p=1e100你有本事給我算啊(°□°)，還不是同樣算不出來）

 

而後，咱們著名的盧卡斯（Lucas）在人羣中站了出來(｀･д･´)說：「讓老子來教你這題」

 

盧卡斯說：

C(n, m) % p  =  C(n / p, m / p) * C(n%p, m%p) % p

對於C(n / p, m / p)，若是n / p 仍是很大，能夠遞歸下去，一直到世界的盡頭

 

 

衆人聞此言，無不驚歎，妙哉！妙哉！

 

（笑死我了o(*≧▽≦)ツ┏━┓拍桌狂笑）

（不要問我證實過程，我不費(´･ω･｀)）

 

 

而後上代碼

1 LL Lucas(LL n, LL m, int p){
2         return m ? Lucas(n/p, m/p, p) * comb(n%p, m%p, p) % p : 1;
3 }

 

 

簡潔易懂<(￣︶￣)>

 

 

 

 

 

 

 

實戰一下吧

 

hdu 5446

http://acm.hdu.edu.cn/showproblem.php?pid=5446

 

題意：

給你三個數n, m, k

第二行是k個數，p1,p2,p3...pk

全部p的值不相同且p都是質數

求C(n, m) % (p1*p2*p3*...*pk)的值

範圍：1≤m≤n≤1e18，1≤k≤10，pi≤1e5，保證p1*p2*p3*...*pk≤1e18

 

 

AC代碼：

 1 #include<cstdio>
 2 typedef long long LL;
 3 const int N = 100000 + 5;
 4 LL mul(LL a, LL b, LL p){//快速乘，計算a*b%p 
 5     LL ret = 0;
 6     while(b){
 7         if(b & 1) ret = (ret + a) % p;
 8         a = (a + a) % p;
 9         b >>= 1;
10     }
11     return ret;
12 }
13 LL fact(int n, LL p){//n的階乘求餘p 
14     LL ret = 1;
15      for (int i = 1; i <= n ; i ++) ret = ret * i % p ;
16     return ret ;
17 }
18 void ex_gcd(LL a, LL b, LL &x, LL &y, LL &d){
19     if (!b) {d = a, x = 1, y = 0;}
20     else{
21         ex_gcd(b, a % b, y, x, d);
22         y -= x * (a / b);
23     }
24 }
25 LL inv(LL t, LL p){//若是不存在，返回-1 
26     LL d, x, y;
27     ex_gcd(t, p, x, y, d);
28     return d == 1 ? (x % p + p) % p : -1;
29 }
30 LL comb(int n, int m, LL p){//C(n, m) % p
31     if (m < 0 || m > n) return 0;
32     return fact(n, p) * inv(fact(m, p), p) % p * inv(fact(n-m, p), p) % p;
33 }
34 LL Lucas(LL n, LL m, int p){
35         return m ? Lucas(n/p, m/p, p) * comb(n%p, m%p, p) % p : 1;
36 }
37 LL china(int n, LL *a, LL *m){//中國剩餘定理 
38     LL M = 1, ret = 0;
39     for(int i = 0; i < n; i ++) M *= m[i];
40     for(int i = 0; i < n; i ++){
41         LL w = M / m[i];
42         //ret = (ret + w * inv(w, m[i]) * a[i]) % M;//這句寫了會WA，用下面那句 
43         ret = (ret + mul(w * inv(w, m[i]), a[i], M)) % M;
44         //由於這裏直接乘會爆long long ,因此我用快速乘(unsigned long long也是爆掉，除非用高精度)
45     }
46     return (ret + M) % M;
47 }
48 int main(){
49     int T, k;
50     LL n, m, p[15], r[15];
51     scanf("%d", &T);
52     while(T--){
53         scanf("%I64d%I64d%d", &n, &m, &k);
54         for(int i = 0; i < k; i ++){
55             scanf("%I64d", &p[i]);
56             r[i] = Lucas(n, m, p[i]);
57         }
58         printf("%I64d\n", china(k, r, p));
59     }
60 }

View Code

 

 

咱們知道題目要求C(n, m) % (p1*p2*p3*...*pk)的值

其實這個就是中國剩餘定理最後算出結果後的最後一步求餘

那C(n, m)至關於之前咱們須要用中國剩餘定理求的值

然而C(n, m)太大，咱們只好先算出

C(n, m) % p1 = r1

C(n, m) % p2 = r2

C(n, m) % p3 = r3

.

.

.

C(n, m) % pk = rk

用Lucas，這些r1,r2,r3...rk能夠算出來

而後又是用中國剩餘定理求答案

 

 

 

 

 

 

 

 

 

全都是套路。。。