组合数取模&&Lucas定理题集
Posted 努力努力再努力x
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了组合数取模&&Lucas定理题集相关的知识，希望对你有一定的参考价值。
题集链接：
https://cn.vjudge.net/contest/231988
解题之前请先了解组合数取模和Lucas定理
A : FZU-2020
输出组合数C(n, m) mod p
(1 <= m <= n <= 10^9, m <= 10^4, m < p < 10^9, p是素数)
由于p较大，不可以打表，直接Lucas求解
1 #include<iostream> 2 using namespace std; 3 typedef long long ll; 4 const int maxn = 1e5 + 10; 5 ll fac[maxn];//阶乘打表 6 void init(ll p)//此处的p应该小于1e5，这样Lucas定理才适用 7 { 8 fac[0] = 1; 9 for(int i = 1; i <= p; i++) 10 fac[i] = fac[i - 1] * i % p; 11 } 12 ll pow(ll a, ll b, ll m) 13 { 14 ll ans = 1; 15 a %= m; 16 while(b) 17 { 18 if(b & 1)ans = (ans % m) * (a % m) % m; 19 b /= 2; 20 a = (a % m) * (a % m) % m; 21 } 22 ans %= m; 23 return ans; 24 } 25 ll inv(ll x, ll p)//x关于p的逆元，p为素数 26 { 27 return pow(x, p - 2, p); 28 } 29 ll C(ll n, ll m, ll p)//组合数C(n, m) % p 30 { 31 if(m > n)return 0; 32 ll up = 1, down = 1;//分子分母; 33 for(int i = n - m + 1; i <= n; i++)up = up * i % p; 34 for(int i = 1; i <= m; i++)down = down * i % p; 35 return up * inv(down, p) % p; 36 } 37 ll Lucas(ll n, ll m, ll p) 38 { 39 if(m == 0)return 1; 40 return C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p; 41 } 42 int main() 43 { 44 ll n, m, p; 45 int T; 46 cin >> T; 47 while(T--) 48 { 49 cin >> n >> m >> p; 50 //init(p); 51 cout<<Lucas(n, m, p)<<endl; 52 } 53 return 0; 54 }
B：HDU-3944
求从顶点到达第n行，第m列的元素的最短路径值（最开始是第0行，第0列）
只能向下或者斜向右。
如图，如果在右侧，可以走红色路线，比如第4行第3列，与此对称的是第4行第1列（从第0行，第0列数起）
这样，我们就可以把左侧的转化成右侧来统一计算。中间的也可以这样计算。
左侧转化成右侧：对于n行m列而言，如果m * 2 < n 那么m = n - m
对于右侧的如何计算呢？
用C(x, y)表示x行y列的值，也就是组合数的值
那么ans = C(n , m) + C(n - 1, m) + C(n - 2, m) + ... + C(m + 1, m) + C(m, m) + C(m - 1, m - 1) + C(m - 2, m - 2) + C(0, 0)
红色部分为m
蓝色部分可以通过组合数公式C(n, m) + C(n, m + 1) = C(n +1, m + 1)来计算
把蓝色部分的C(m, m)写成C(m + 1, m + 1)
那么最后两项就可以通过上面那个公式合并成C(m + 2, m + 1)
不断合并，最后得到C(n, m) + C(n, m + 1) = C(n +1, m + 1)
所以答案就是C(n + 1, m + 1) + m % p;
由于计算量大，p的范围是10000内的素数，而需要计算10^6个组合数，所以需要先打表，把10000内的素数筛选出来，打表每个阶乘模上素数的值，还有逆元，这样计算组合数才不会超时，而且数组不能开long long，会超时
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef long long ll; 4 const int maxn = 1e4 + 10; 5 int fac[maxn][maxn], inv[maxn][maxn];//阶乘打表 逆元打表 6 int prime[maxn], pth[maxn]; 7 int pow(int a, int b, int m) 8 { 9 int ans = 1; 10 a %= m; 11 while(b) 12 { 13 if(b & 1)ans = (ans % m) * (a % m) % m; 14 b /= 2; 15 a = (a % m) * (a % m) % m; 16 } 17 ans %= m; 18 return ans; 19 } 20 bool is_prime[maxn]; 21 void init() 22 { 23 int n = 10000, cnt = 0; 24 for(int i = 2; i <= n; i++)is_prime[i] = 1; 25 for(int i = 2; i <= n; i++) 26 { 27 if(is_prime[i]) 28 { 29 prime[++cnt] = i; 30 pth[i] = cnt; 31 for(int j = i * i; j <= n; j += i)is_prime[j] = 0; 32 } 33 } 34 //cout<<cnt<<endl; 35 for(int i = 1; i <= cnt; i++) 36 { 37 fac[i][0] = inv[i][0] = 1; 38 for(int j = 1; j < prime[i]; j++)//大于i的阶乘模上i均为0 39 { 40 fac[i][j] = fac[i][j - 1] * j % prime[i]; 41 inv[i][j] = pow(fac[i][j], prime[i] - 2, prime[i]); 42 } 43 } 44 } 45 46 int C(ll n, int m, int p)//组合数C(n, m) % p 47 { 48 if(m > n)return 0; 49 if(m == n)return 1; 50 int t = pth[p]; 51 return fac[t][n] * (inv[t][m] * inv[t][n - m] % p) % p; 52 } 53 54 int Lucas(int n, int m, int p) 55 { 56 //cout<<n<<" "<<m<<" "<<p<<endl; 57 if(m == 0)return 1; 58 return C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p; 59 } 60 61 int main() 62 { 63 init(); 64 int n, m, p; 65 int cases = 0; 66 while(scanf("%d%d%d", &n, &m, &p) != EOF) 67 { 68 if(m <= n / 2) 69 { 70 m = n - m; 71 } 72 ll ans = (Lucas(n + 1, m + 1, p) + m) % p; 73 printf("Case #%d: %d\\n", ++cases, ans); 74 } 75 return 0; 76 }
C:ZOJ-3557
在n个元素的集合中取出lm和不相邻元素的种数
插空法，由于不相邻，那就把这m个先取出来，放在剩下的n-m个数字的两侧的空位上，那么这m个一定不相邻。
等价于n-m+1中放m个元素
C（n - m + 1, m） % p
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 inv(ll x, ll p)//x关于p的逆元，p为素数 20 { 21 return pow(x, p - 2, p); 22 } 23 ll C(ll n, ll m, ll p)//组合数C(n, m) % p 24 { 25 if(m > n)return 0; 26 ll up = 1, down = 1;//分子分母; 27 for(int i = n - m + 1; i <= n; i++)up = up * i % p; 28 for(int i = 1; i <= m; i++)down = down * i % p; 29 return up * inv(down, 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 } 36 int main() 37 { 38 ll n, m, p; 39 while(cin >> n >> m >> p) 40 { 41 if(n - m + 1 < m) 42 cout<<"0"<<endl; 43 else cout<<Lucas(n - m + 1, m, p)<<endl; 44 } 45 return 0; 46 }
D在F介绍完之后再介绍
E：hdu-4349
求C(n,0),C (n,1),C (n,2)...C (n,n)中奇数的数目
首先：组合数奇偶性判断公式：n & m == m
当然本题要判断的组合数很多，所以不能用上述结论，只能另辟蹊径。由于是判断奇偶性，那么就是判断 
是否为1，利用Lucas定理，先把
和
化为二进制，这样它们都是01序列了。我们又知道 
。这样中为0的地方对应的中的位置只有一种可能，那就是0。
这样我们可以不用管中为0的地方，只考虑中为1的位置，可以看出，中为1的位置对应的中为0或1，其结果都是1，这样答案就是：1<<(二进制表示中1的个数)
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 int lowbit(int x) 7 { 8 return x & (-x); 9 } 10 int main() 11 { 12 ll n, m, p; 13 while(cin >> n) 14 { 15 int tot = 0; 16 while(n) 17 n -= lowbit(n), tot++; 18 cout<<(1<<tot)<<endl; 19 } 20 return 0; 21 }
F：Gym - 100633J
扩展Lucas定理模板
这里的p不是素数，用中国剩余定理合并同余方程组。
详细介绍请看全文最开始的传送门
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 }
D：HDU-4373
个for循环嵌套，有两种形式，第一类从1开始到，第二类从上一层循环当前数开始到，第一层一定是第一种类型，求总的循环的次数对364875103取余的结果。
首先可以看出，每一个第一类循环都是一个新的开始，与前面的状态无关，所以可以把个嵌套分为几个不
同的部分，每一个部分由第一类循环开始，最终结果相乘即可。剩下的就是第二类循环的问题，假设一个
层循环，最大到，分析一下得到如下结果
（1）只有一层，则循环次数为
（2）只有前两层，则循环次数为
（3）只有前三层，则循环次数为
由此得到结论：第的循环次数为
由于364875103=97 * 3761599
用中国剩余定理合并，但是不能直接套上述模板，会超时，上述模板是对于不同的合数p分解成质数的次方，算起来很麻烦，此处364875103就等于两个质数的乘积
这里a1 = 97, a2 = 3761599所以M1 = 3761599， M2 = 97， 求出M1关于m1的逆元，M2关于m2的逆元
然后之需要求出组合数mod97的值为a1，组合数mod3761599的值为a2
ans按照上述公式就可以直接算出来了。
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef long long ll; 4 const int maxn = 1e6 + 10; 5 ll pow(ll a, ll b, ll m) 6 { 7 ll ans = 1; 8 a %= m; 9 while(b) 10 { 11 if(b & 1)ans = (ans % m) * (a % m) % m; 12 b /= 2; 13 a = (a % m) * (a % m) % m; 14 } 15 ans %= m; 16 return ans; 17 } 18 const ll mod1 = 97; 19 const ll mod2 = 3761599; 20 const ll mod = mod1 * mod2; 21 ll fac1[mod1 + 10], fac2[mod2 + 10]; 22 ll inv1, inv2; 23 void init() 24 { 25 fac1[0] = fac2[0] = 1; 26 for(int i = 1; i < mod1; i++)fac1[i] = fac1[i - 1] * i % mod1; 27 for(int i = 1; i < mod2; i++)fac2[i] = fac2[i - 1] * i % mod2; 28 inv1 = pow(mod2, mod1 - 2, mod1);//mod1 模上mod2的逆元 29 inv2 = pow(mod1, mod2 - 2, mod2);//mod2 模上mod1的逆元 求出在之后的中国剩余定理合并的时候不用计算 30 } 31 32 ll C(ll n, ll m, ll p, ll fac[]) 33 { 34 if(m > n)return 0; 35 return fac[n] * pow(fac[m] * fac[n - m], p - 2, p) % p; 36 } 37 38 ll Lucas(ll n, ll m, ll p, ll fac[]) 39 { 40 if(m == 0)return 1; 41 return C(n % p, m % p, p, fac) * Lucas(n / p, m / p, p, fac) % p; 42 } 43 int main() 44 { 45 ll n, m, k, cases = 0; 46 ll a[20]; 47 int T; 48 Lucas定理及组合数取模