在数论中我们经常遇到求组合数的问题,比如有n个球放入m个箱子里(m > n),每个箱子只能放一个球,问有多少种方案?类似这种题目就需要用到组合数来求解
学习该内容需学会之前讲的逆元和中国剩余定理。跳转链接:逆元 中国剩余定理
那么求组合数都有什么方法呢?
一.其实早在古代,伟大的先人就解决了求解组合数的问题,那就是伟大的杨辉三角
我们发现
1
1,1
1,2,1
1,3,3,1
1,4,6,4,1
1,5,10,10,5,1
有没有发现,C(n,m)的值就是第n行m列的值(下标从0开始)
虽然这样可以求组合数,但是它的复杂度是O(n^2)太大了,只能求解1e3以内的组合数
二.我们可以进一步优化
由于已知C(n,m)= n!/ (n - m)! ✖️m!
我们只需要预处理出前n的阶乘与前n阶乘的逆元,每次就可以O(1)的时间求出组合数的值
因此,我们只需要耗费O(n)的复杂度预处理就可以了
因此它可以求1e6以内的组合数
三.卢卡斯定理
卢卡斯定理:C(n,m)% p = C(n / p ,m / p)✖️C(n % p ,m % p) % p
当n和m特别大时题目中必须有取余操作,要不会爆longlong,若p小于1e7时卢卡斯定理就派上了用场。
附上代码:
ll fact[N + 5];
ll inv[N + 5];
//省略中间的扩展欧几里得求逆元
void init () {
fact[0] = fact[1] = 1;
for (int i = 2; i <= N; i ++) {
fact[i] = fact[i - 1] * i % mod;
}
inv[N] = mod_reverse(fact[N], mod); //调用扩展欧几里得求逆元
for (int i = N -1; i >= 0; i--) {
inv[i] = inv[i + 1] * (i + 1) % mod;
}
return;
}
ll C (ll n, ll m) {
if (n < m || m < 0) return 0;
if (n < mod) return fact[n] * inv[m] % mod * inv[n - m] % mod;
return C (n / mod, m / mod) * C (n % mod, m % mod) % mod;
}
以上代码只适用于p<1e7且p为素数的情况,若p为合数呢?
此时就需要用到扩展卢卡斯定理了(此前需学会中国剩余定理)
四.扩展卢卡斯定理:
由于p为合数,因此我们可以将p分解为p1a1✖️p2a2✖️…pk^ak的形式
因此我们可以利用卢卡斯定理求得:
ans ≡ x1 (mod p1^a1)
ans ≡ x2 (mod p2^a2)
…
ans ≡ xk (mod pk^ak)
不同质数的幂次方之间肯定两两互质,通过求得每组的解xi,再利用中国剩余定理就可以合并得到最终答案
每组解xi的求法:C(n,m)% (pi^ai)
因此我们需要获得剔除pi因子后的乘积n!,m!的逆元,(n - m)!的逆元的乘积去取模(pi^ai)并加入被剔除的pi的贡献即可
那么我们如何求n!剔除pi因子的乘积呢?
例如19! % (3^2) = (1✖️2✖️4✖️5✖️7✖️8)✖️(10✖️11✖️13✖️14✖️16✖️17)✖️(19)✖️(1✖️2✖️3✖️4✖️5✖️6)
()内的代表一个分组,每个分组在取余pi^ai后结果相同,因此只需暴力一个分组,利用快速幂求解即可,对于(19)这种不足一个分组的,可以暴力单度计算。对于(1✖️2✖️3✖️4✖️5✖️6),可提前预处理出来
最后xi = (n)! % (pi^ai)✖️inv((m)! % (pi^ai), pi^ai)✖️inv((m)! % (pi^ai), piai)✖️pik即可(其中k为pi在C(n, m )中出现的次数)
最后利用中国剩余定理将xi合并即可
附上代码:
ll ksm(ll a, ll b, ll mod) {} //省略快速幂
ll extend_gcd(ll a, ll b, ll &x, ll &y) {} //省略扩展欧几里得
ll mod_reverse (ll a, ll n) {} //省略逆元
//用于求n!%(pi^ai) 其中pk =(pi^ai)
ll mul(ll n,ll pi,ll pk) {
if(!n) return 1;
ll ans=1;
if(n / pk) {
for(ll i = 2;i < pk; i++)
if(i % pi) ans = ans * i % pk;
ans = ksm(ans, n / pk, pk);
}
for(ll i = 2;i <= n % pk; i++)
if(i % pi) ans = ans * i % pk;
return ans * mul(n / pi, pi, pk) % pk;
}
//求xi
ll C(ll n, ll m, ll mod, ll pi, ll pk) {
if(m > n) return 0;
ll a = mul(n, pi, pk),b = mul(m, pi, pk),c = mul(n-m, pi, pk);
ll k=0, ans;
for(ll i = n; i; i /= pi) k += i / pi;
for(ll i = m; i; i /= pi) k -= i / pi;
for(ll i = n - m; i; i /= pi) k -= i / pi;
ans = a * mod_reverse(b, pk) % pk * mod_reverse(c, pk) % pk * ksm(pi, k, pk) % pk;
return ans * (mod / pk) % mod * mod_reverse(mod / pk, pk) % mod;
}
//扩展卢卡斯
ll exLucas(ll n, ll m, ll mod) {
ll ans = 0;
for(ll x = mod, i = 2; i <= mod; i++)
if(x % i == 0) {
ll pk = 1;
while(x % i == 0) {
pk *= i;
x /= i;
}
ans = (ans + C(n, m, mod, i, pk)) % mod;
}
return ans;
}
完工~
如果有写的不对或者不全面的地方 可通过主页的联系方式进行指正,谢谢