组合数取模方法总结(Lucas定理介绍)

1.当n,m都很小的时候可以利用杨辉三角直接求。
C(n,m)=C(n-1,m)+C(n-1,m-1);

 

2、n和m较大,但是p为素数的时候

Lucas定理是用来求 c(n,m) mod p,p为素数的值

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

 

也就是Lucas(n,m)%p=Lucas(n/p,m/p)*C(n%p,m%p)%p

 

 

求上式的时候,Lucas递归出口为m=0时返回1

求C(n%p, m%p)%p的时候,此处写成C(n, m)%p(p是素数,n和m均小于p)

C(n, m)%p = n! / (m ! * (n - m )!) % p = n! * mod_inverse[m! * (n - m)!, p] % p

由于p是素数,有费马小定理可知,m! * (n - m)! 关于p的逆元就是m! * (n - m)!的p-2次方。

 

p较小的时候预处理出1-p内所有阶乘%p的值,然后用快速幂求出逆元,就可以求出解。p较大的时候只能逐项求出分母和分子模上p的值,然后通过快速幂求逆元求解。

P较大,不打表:

 

 1 ll pow(ll a, ll b, ll m)
 2 {
 3     ll ans = 1;
 4     a %= m;
 5     while(b)
 6     {
 7         if(b & 1)ans = (ans % m) * (a % m) % m;
 8         b /= 2;
 9         a = (a % m) * (a % m) % m;
10     }
11     ans %= m;
12     return ans;
13 }
14 ll inv(ll x, ll p)//x关于p的逆元,p为素数
15 {
16     return pow(x, p - 2, p);
17 }
18 ll C(ll n, ll m, ll p)//组合数C(n, m) % p
19 {
20     if(m > n)return 0;
21     ll up = 1, down = 1;//分子分母;
22     for(int i = n - m + 1; i <= n; i++)up = up * i % p;
23     for(int i = 1; i <= m; i++)down = down * i % p;
24     return up * inv(down, p) % p;
25 }
26 ll Lucas(ll n, ll m, ll p)
27 {
28     if(m == 0)return 1;
29     return C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p;
30 }

 

P较小,打表:

 1 const int maxn = 1e5 + 10;
 2 ll fac[maxn];//阶乘打表
 3 void init(ll p)//此处的p应该小于1e5,这样Lucas定理才适用
 4 {
 5     fac[0] = 1;
 6     for(int i = 1; i <= p; i++)
 7         fac[i] = fac[i - 1] * i % p;
 8 }
 9 ll pow(ll a, ll b, ll m)
10 {
11     ll ans = 1;
12     a %= m;
13     while(b)
14     {
15         if(b & 1)ans = (ans % m) * (a % m) % m;
16         b /= 2;
17         a = (a % m) * (a % m) % m;
18     }
19     ans %= m;
20     return ans;
21 }
22 ll inv(ll x, ll p)//x关于p的逆元,p为素数
23 {
24     return pow(x, p - 2, p);
25 }
26 ll C(ll n, ll m, ll p)//组合数C(n, m) % p
27 {
28     if(m > n)return 0;
29     return fac[n] * inv(fac[m] * fac[n - m], p) % p;
30 }
31 ll Lucas(ll n, ll m, ll p)
32 {
33     if(m == 0)return 1;
34     return C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p;
35 }

 

3、n,m较大且p不为素数的时候

扩展Lucas定理:

 

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 typedef long long ll;
 4 const int maxn = 1e6 + 10;
 5 const int mod = 1e9 + 7;
 6 ll pow(ll a, ll b, ll m)
 7 {
 8     ll ans = 1;
 9     a %= m;
10     while(b)
11     {
12         if(b & 1)ans = (ans % m) * (a % m) % m;
13         b /= 2;
14         a = (a % m) * (a % m) % m;
15     }
16     ans %= m;
17     return ans;
18 }
19 ll extgcd(ll a, ll b, ll& x, ll& y)
20 //求解ax+by=gcd(a, b)
21 //返回值为gcd(a, b)
22 {
23     ll d = a;
24     if(b)
25     {
26         d = extgcd(b, a % b, y, x);
27         y -= (a / b) * x;
28     }
29     else x = 1, y = 0;
30     return d;
31 }
32 ll mod_inverse(ll a, ll m)
33 //求解a关于模上m的逆元
34 //返回-1表示逆元不存在
35 {
36     ll x, y;
37     ll d = extgcd(a, m, x, y);
38     return d == 1 ? (m + x % m) % m : -1;
39 }
40 
41 ll Mul(ll n, ll pi, ll pk)//计算n! mod pk的部分值  pk为pi的ki次方
42 //算出的答案不包括pi的幂的那一部分
43 {
44     if(!n)return 1;
45     ll ans = 1;
46     if(n / pk)
47     {
48         for(ll i = 2; i <= pk; i++) //求出循环节乘积
49             if(i % pi)ans = ans * i % pk;
50         ans = pow(ans, n / pk, pk); //循环节次数为n / pk
51     }
52     for(ll i = 2; i <= n % pk; i++)
53         if(i % pi)ans = ans * i % pk;
54     return ans * Mul(n / pi, pi, pk) % pk;//递归求解
55 }
56 
57 ll C(ll n, ll m, ll p, ll pi, ll pk)//计算组合数C(n, m) mod pk的值 pk为pi的ki次方
58 {
59     if(m > n)return 0;
60     ll a = Mul(n, pi, pk), b = Mul(m, pi, pk), c = Mul(n - m, pi, pk);
61     ll k = 0, ans;//k为pi的幂值
62     for(ll i = n; i; i /= pi)k += i / pi;
63     for(ll i = m; i; i /= pi)k -= i / pi;
64     for(ll i = n - m; i; i /= pi)k -= i / pi;
65     ans = a * mod_inverse(b, pk) % pk * mod_inverse(c, pk) % pk * pow(pi, k, pk) % pk;//ans就是n! mod pk的值
66     ans = ans * (p / pk) % p * mod_inverse(p / pk, pk) % p;//此时用剩余定理合并解
67     return ans;
68 }
69 
70 ll Lucas(ll n, ll m, ll p)
71 {
72     ll x = p;
73     ll ans = 0;
74     for(ll i = 2; i <= p; i++)
75     {
76         if(x % i == 0)
77         {
78             ll pk = 1;
79             while(x % i == 0)pk *= i, x /= i;
80             ans = (ans + C(n, m, p, i, pk)) % p;
81         }
82     }
83     return ans;
84 }
85 
86 int main()
87 {
88     ll n, m, p;
89     while(cin >> n >> m >> p)
90     {
91         cout<<Lucas(n, m, p)<<endl;
92     }
93     return 0;
94 }

 

越努力,越幸运
posted @ 2018-05-28 16:01  _努力努力再努力x  阅读(8056)  评论(7编辑  收藏  举报