（数论十）卢卡斯定理与扩展卢卡斯

​ 在数论中我们经常遇到求组合数的问题，比如有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分解为p1^a1✖️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), pi^ai)✖️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;
}

​ 完工～

转载请注明出处！！！

如果有写的不对或者不全面的地方 可通过主页的联系方式进行指正，谢谢