组合数取模&&Lucas定理题集

题集链接：

https://cn.vjudge.net/contest/231988

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
", ++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     init();
49     cin >> T;
50     while(T--)
51     {
52         cin >> n >> m;
53         cin >> k;
54         for(int i = 0; i < k; i++)cin >> a[i];
55         a[k] = m;
56         ll ans = 1;
57         for(int i = 1; i <= k; i++)
58         {
59             ll t = a[i] - a[i - 1];
60             ll a1 = Lucas(n + t - 1, t, mod1, fac1);
61             ll a2 = Lucas(n + t - 1, t, mod2, fac2);
62             ll a3 = (a1 * mod2 % mod * inv1 + a2 * mod1 % mod * inv2) % mod;//中国剩余定理合并同余方程组
63             ans = ans * a3 % mod;
64         }
65         cout<<"Case #"<<++cases<<": "<<ans<<endl;
66     }
67     return 0;
68 }