这是一个模板。卢卡斯定理用于大组合数求模运算。

        首先引入卢卡斯定理。

【卢卡斯定理】对于整数$n，m（m>0）和质数p，$若$n=sp+q , m=tp+r$，则有$C^{n}_m mod p=C^s_tC^q_r mod p$。

        证明略。该定理可以用于快速求大组合数对某个质数的模。这里有模板题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#include<bits/stdc++.h>

using namespace std;
long long n, m, p, op[100005];

long long qPow(long long a, long long b, long long p) {//快速幂
    long long  s = a, ans = 1ll;
    while (b) {
        if (b & 1)ans = ans * s % p;
        s = s * s % p, b >>= 1;
    }
    return ans;
}

long long C(long long a, long long b, long long p) {
    if (b < a)return 0ll;//仍然注意边界处理
    long long s = op[a] * op[b - a] % p;//根据定义式，求分子
    return op[b] * qPow(s, p - 2, p) % p;//乘法逆元处理
}

long long lucas(long long x, long long y, long long p) {
    if (x == 0)return 1ll;//注意处理边界
    return lucas(x / p, y / p, p) * C(x % p, y % p, p) % p;//后者直接求，前者递归求
}

int main() {
    int t;
    cin >> t;
    while (t--) {
        cin >> n >> m >> p;
        op[0] = op[1] = 1ll;
        for (int i = 2; i <= p; i++)op[i] = op[i - 1] * i % p;//预处理阶乘
        cout << lucas(m, m + n, p) << endl;
    }
    return 0;
}

扩展卢卡斯

        如果模数不是质数，卢卡斯定理便不再适用，这个时候如何求解组合数的模？需要用到扩展卢卡斯，只不过这一个算法与卢卡斯定理并没有什么关系，它需要前缀知识：扩展欧几里得算法和中国剩余定理（CRT）。
        对于模数mod，根据整数分解唯一性定理，可以将其分解成下面的形式：

        如果可以求出$C_n^m$对$p^s$的模，那么可以利用CRT合并答案，问题转化为如何求组合数对质数的幂次的模。
        由于分母不一定与质数幂次互质，其逆元不一定存在，无法用逆元来求模。这里可以根据组合数定义：

        将分子分母的所有因子p单独提出，这样就可以用逆元了，最后再乘上这些p因子乘积对$p^s$的模即可。下面考虑如何在剔除因子p的情况下快速求解阶乘。
        现在假设要求解$x!$，质数为$p$，模数为$p^s$。将x分成几块：$1$到$p^s$是第一块，$p^s+1$到$2p^s$是第二块，依次类推，最后剩余的数暂且不去管它。可以发现，每一块都含有相同数量的p的倍数，并且对于每一块，将这些p的倍数剔除，剩下的数乘起来模$p^s$的值都是相同的。对于x来说，这样的块有$\lfloor \frac {x} {p^s} \rfloor$块，这样只需求解第一块除去p的倍数的数的乘积，再用一个快速幂即可。对于剩下的数，将它们暴力乘起来。
        现在已经处理了所有非p的倍数，对于p的倍数，将它们依次除去因子p（前文已经说明需要剔除因子p），可以得到一个阶乘，容易知道这个阶乘是$\lfloor \frac {x} {p}\rfloor !$，递归求解即可，递归终止条件是$0!=1$。

1
2
3
4
5
6
7
8
ll fac(ll a, ll p, ll pk) {//求阶乘：a!，p是质数，pk是p的幂次，这里定义long long 为ll
    ll ans = 1ll;
    if (a == 0)return 1ll;
    for (ll i = 2; i <= pk; i++)if (i % p)ans = ans * i % pk;//第一块
    ans = qPow(ans, a / pk, pk);//qPow是快速幂
    for (ll i = 2; i <= a % pk; i++)if (i % p)ans = ans * i % pk;//暴力处理剩余的
    return ans * fac(a / p, p, pk) % pk;//递归
}

        那现在如何知道$x!$中有多少p的因子呢？首先p的倍数有$\lfloor \frac {x} {p}\rfloor$个，将它们全部除以p，得到一个阶乘$\lfloor \frac {x} {p}\rfloor !$，再计入它的p因子数目即可。若令$f(x)$表示$x!$中因子p的数目，由此得到递推式：

        这样组合数求解函数就可以写出了：

1
2
3
4
5
6
7
8
ll C(ll a, ll b, ll p, ll pk) {
    ll x = fac(a, p, pk), y = fac(b, p, pk), z = fac(a - b, p, pk);//先求阶乘
    ll k = 0;
    for (ll i = a; i; i /= p)k += i / p;//递归求p因子数目
    for (ll i = b; i; i /= p)k -= i / p;
    for (ll i = a - b; i; i /= p)k -= i / p;//这是k就是提出来的p因子数目
    return x * inv(y, pk) % pk * inv(z, pk) % pk * qPow(p, k, pk) % pk;//inv是逆元
}

        再来一步中国剩余定理即可：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
ll CRT(ll a, ll m, ll mod) {//中国剩余定理
    return a * inv(mod / m, m) % mod * (mod / m) % mod;
}

ll exLucas(ll n, ll m, ll mod) {//扩展卢卡斯主函数
    ll ans = 0ll, tmp = mod, pk;
    for (ll i = 2; i * i <= tmp; i++) {
        if (tmp % i == 0) {
            pk = 1ll;
            while (tmp % i == 0)pk = pk * i, tmp /= i;
            ans = (ans + CRT(C(n, m, i, pk), pk, mod)) % mod;
        }
    }
    if (tmp > 1)ans = (ans + CRT(C(n, m, tmp, tmp), tmp, mod)) % mod;
    return ans;
}

        由于模数不一定是质数，故采用扩展欧几里得的方式来求逆元：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
ll exGcd(ll a, ll b, ll &x, ll &y) {//扩展欧几里得
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    ll r = exGcd(b, a % b, y, x);
    y -= a / b * x;
    return r;
}

ll inv(ll a, ll b) {//求逆元
    ll x, y;
    exGcd(a, b, x, y);
    return (x % b + b) % b;//保证非负
}

        附模板题全代码：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#include<bits/stdc++.h>

using namespace std;
typedef long long ll;

ll exGcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    ll r = exGcd(b, a % b, y, x);
    y -= a / b * x;
    return r;
}

ll qPow(ll a, ll b, ll mod) {
    ll ans = 1ll, sta = a;
    while (b) {
        if (b & 1)ans = ans * sta % mod;
        sta = sta * sta % mod, b >>= 1;
    }
    return ans;
}

ll inv(ll a, ll b) {
    ll x, y;
    exGcd(a, b, x, y);
    return (x % b + b) % b;
}

ll CRT(ll a, ll m, ll mod) {
    return a * inv(mod / m, m) % mod * (mod / m) % mod;
}

ll fac(ll a, ll p, ll pk) {
    ll ans = 1ll;
    if (a == 0)return 1ll;
    for (ll i = 2; i <= pk; i++)if (i % p)ans = ans * i % pk;
    ans = qPow(ans, a / pk, pk);
    for (ll i = 2; i <= a % pk; i++)if (i % p)ans = ans * i % pk;
    return ans * fac(a / p, p, pk) % pk;
}

ll C(ll a, ll b, ll p, ll pk) {
    ll x = fac(a, p, pk), y = fac(b, p, pk), z = fac(a - b, p, pk);
    ll k = 0;
    for (ll i = a; i; i /= p)k += i / p;
    for (ll i = b; i; i /= p)k -= i / p;
    for (ll i = a - b; i; i /= p)k -= i / p;
    return x * inv(y, pk) % pk * inv(z, pk) % pk * qPow(p, k, pk) % pk;
}

ll exLucas(ll n, ll m, ll mod) {
    ll ans = 0ll, tmp = mod, pk;
    for (ll i = 2; i * i <= tmp; i++) {
        if (tmp % i == 0) {
            pk = 1ll;
            while (tmp % i == 0)pk = pk * i, tmp /= i;
            ans = (ans + CRT(C(n, m, i, pk), pk, mod)) % mod;
        }
    }
    if (tmp > 1)ans = (ans + CRT(C(n, m, tmp, tmp), tmp, mod)) % mod;
    return ans;
}

int main() {
    ll n, m, mod;
    cin >> n >> m >> mod;
    cout << exLucas(n, m, mod);
    return 0;
}