数论 · exLucas 定理

pldzy · 2021-12-24 21:30:24 · 个人记录

前言

终于看懂了！！！连夜丢掉卡特兰来做了两道 exLucas 的题目。

定理

作为一个没什么用的铺垫：数论 · Lucas 定理

求解 C_n^m \%p（不保证 p 为质数）。

证明

1

先把 p 拆分一下：p = p_1^{a_1} * p_2^{a_2} * \cdots p_k^{a_k}。

然后逐一求解 C_n^m \% p_i^{a_i}，

最后用 CRT 求解 \begin{cases}x\equiv{b_1}\pmod{p_1^{a_1}}\\x\equiv{b_2}\pmod{p_2^{a_2}}\\\cdots\\x\equiv{b_k}\pmod{p_k^{a_k}}\end{cases} 的最小正整数解 x 即可。

2

假设求解 C_n^m \% p_i^{k}，设 p_k=p_i^k。

有 C_n^m=\dfrac{n!}{m!(n - m)!}，

所以就变成求解 \dfrac{n!\%p_k}{m!\%p_k(n - m)!\%p_k}。

但是因为分母部分需要逆元，有可能 m!\%p_k 或 (n - m)!\%p_k 与 p_k 不互质，所以我们要先除去他们中所有的 p_i，然后取模完之后再乘回去即可。

假设 x 为 n! 中包含的 p_i 的个数，y、z 同理。

则变为求解 (\dfrac{\dfrac{n!}{p_i^x} \% p_k}{\dfrac{m!}{p_i^y} \% p_k\dfrac{(n-m)!}{p_i^z} \% p_k} * p_i^{x-y-z}) \% p_k

3

考虑求解 \dfrac{n!}{p_i^x}\% p_k。

先举个例子：

设 n=19,\ p_i = 3,\ k = 2。

先把 n! 中 p_i 的倍数先提取出来再合并一下：

按照 2 的后半部分所述，3^6 可以先忽略掉。

而柿子的后半部分 (1* 2* 3* 4* 5* 6)，也就是 \left\lfloor\dfrac{n}{p_i}\right\rfloor! 我们可以用递归再次求解。

对于柿子的前半部分，会发现它对 p_i^k 取余会有一个周期：

所以只用求出第一个周期：res=\prod_{i=1}^{p_k} i\ ( p_i \nmid i)，然后再用快速幂求出 res^{\left\lfloor\frac{n}{p_k}\right\rfloor} 即可。

然后就会发先剩下个 19，这些数直接暴力乘进 res（进行完快速幂之后）即可。

4

关于 2 中 x,\ y,\ z 的计算。

小奥中学过试除法（忘了叫啥），可以求解类似“100 的阶乘中有多少个因子 5”的问题。

这里直接除然后求出 x-y-z 即可。

代码如下：

int al = n;
while(al)
    k += al / pi, al /= pi;
al = m;
while(al)
    k -= al / pi, al /= pi;
al = n - m;
while(al)
    k -= al / pi, al /= pi;

最后用 扩展欧几里得 求出 \dfrac{m!}{p_i^y} \% p_k\dfrac{(n-m)!}{p_i^z} \% p_k 的逆元即可。

5

关于使用 CRT 合并。

每次求完 C_n^m \% p_i^{k}=b_i，就有了 x\equiv{b_i}\pmod{p_i^{a_i}}。

先给代码：ans += bi * inv(p / pi, pi) % p * (p / pi) % p

bi * inv(p / pi, pi) % p 是为了让它在取模 p_i 的情况下余数为 b_i；

* (p / pi) % p 是为了让它在取模其他质数的情况下余数为 0。

详见 CRT。

代码

#include<algorithm>
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;

#define int long long
#define rep(i, a, b) for(int i = a; i <= b; ++i)

int n, m, p;

inline int power(int x, int t, int mod)
{
    int res = 1;
    while(t >= 1)
    {
        if(t & 1)
            res = res * x % mod;
        t /= 2, x = x * x % mod;
    }
    return res;
}

inline int fac(int nw, int pi, int pk)
{
    if(!nw) return 1;
    int res = 1;
    rep(i, 2, pk)
        if(i % pi) res = res * i % pk;
    res = power(res, nw / pk, pk);
    rep(i, 2, nw % pk)
        if(i % pi) res = res * i % pk;
    return res * fac(nw / pi, pi, pk) % pk;
}

int x, y;

inline void exGcd(int a, int b)
{
    if(!b)
    {
        x = 1, y = 0;
        return;
    }
    exGcd(b, a % b);
    int tmp = x;
    x = y, y = tmp - a / b * y;
    return;
}

inline int inv(int nw, int mod)
{
    exGcd(nw, mod);
    return (x + mod) > mod ? x : x + mod;
}

inline int C(int pi, int pk)
{
    int fn = fac(n, pi, pk), fm = fac(m, pi, pk), fnm = fac(n - m, pi, pk);
    int k = 0;
    int al = n;
    while(al)
        k += al / pi, al /= pi;
    al = m;
    while(al)
        k -= al / pi, al /= pi;
    al = n - m;
    while(al)
        k -= al / pi, al /= pi;
    return fn * inv(fm, pk) % pk * inv(fnm, pk) % pk * power(pi, k, pk) % pk;
}

inline int Crt(int bi, int mi)
{
    return bi * inv(p / mi, mi) % p * (p / mi) % p;
}

inline int exLucas()
{
    int pk, res = 0;
    int lim = ceil(sqrt(p)), tmp = p;
    rep(i, 2, lim)
    {
        if(tmp % i) continue;
        pk = 1;
        while(!(tmp % i)) 
            tmp /= i, pk *= i;
        res = (res + Crt(C(i, pk), pk)) % p;
    } 
    if(tmp > 1)
        res = (res + Crt(C(tmp, tmp), tmp)) % p;
    return res;
}

signed main()
{
    scanf("%lld %lld %lld", &n, &m, &p);
    printf("%lld\n", exLucas());
    return 0;
}

——End——