阶乘取模：从入门到赛场实战

我们今天要解决的核心问题是计算 $n! \pmod{p}$ 的值。其中 $n$ 和 $p$ 可能都很大。

为什么这个问题重要？因为组合计数问题（比如计算 $\binom{n}{k} \pmod p$）的核心就是阶乘及其逆元的计算。搞不定阶乘取模，组合计数基本就没法做。

零、 前置准备：快速幂与模逆元

在深入阶乘之前，确保你已经熟练掌握了两个基本工具：

1. 快速幂 (Exponentiation by Squaring): 高效计算 $a^b \pmod{p}$。
2. 模逆元 (Modular Inverse): 找到一个整数 $x$ 使得 $ax \equiv 1 \pmod{p}$。

如果 $p$ 是素数，根据费马小定理，$a^{p-2} \pmod p$ 就是 $a$ 模 $p$ 的逆元（前提是 $a$ 不被 $p$ 整除）。如果 $p$ 不一定是素数，但 $\gcd(a, p) = 1$，则需要用扩展欧几里得算法 (Extended Euclidean Algorithm, exgcd) 求解 $ax + py = 1$ 中的 $x$。

一、 简单情况：p 为素数，且 n < p

这是最简单的情况。因为 $n < p$ 且 $p$ 是素数，所以 $1, 2, \ldots, n$ 中没有任何一个是 $p$ 的倍数。

直接计算即可：

$$n! \pmod p = (1 \times 2 \times \dots \times n) \pmod p $$

实现非常直接：

#include <bits/stdc++.h>

using namespace std;
using ll = long long;

// 快速幂 a^b mod p
ll power(ll a, ll b, ll p) {
    ll res = 1;
    a %= p;
    while (b > 0) {
        if (b & 1) res = (res * a) % p;
        a = (a * a) % p;
        b >>= 1;
    }
    return res;
}

// 求解 n! mod p, 当 n < p 且 p 是素数
ll fact_mod_simple(int n, ll p) {
    // 输入：n, p (p为素数, n < p)
    // 输出：n! mod p
    if (n >= p) return 0; // 理论上不会进入，因为前提是 n < p
    ll res = 1;
    for (int i = 1; i <= n; ++i) {
        res = (res * i) % p;
    }
    return res;
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);

    int n;
    ll p;
    // 题意：计算 n! mod p，保证 p 是素数且 n < p
    cin >> n >> p; 
    cout << fact_mod_simple(n, p) << endl;

    return 0;
}
// 时间复杂度: O(n)
// 空间复杂度: O(1)

这个场景虽然简单，但它是后面更复杂情况的基础。

二、 理论基石：威尔逊定理 (Wilson's Theorem)

在讨论 $n \ge p$ 的情况前，我们必须了解威尔逊定理。

威尔逊定理: 对于一个素数 $p$，有 $(p-1)! \equiv -1 \pmod p$。反之，若 $(n-1)! \equiv -1 \pmod n$ 且 $n > 1$，则 $n$ 是素数。

这个定理本身不直接计算 $n! \pmod p$ (除非 $n = p-1$)，但它揭示了阶乘与素数模运算之间的深刻联系。它在某些特定问题的证明或推导中非常有用，例如后面我们会看到它在计算 $n! / p^{v_p(n!)} \pmod p$ 时的一个关键应用。

三、 核心挑战：p 为素数，且 n >= p

当 $n \ge p$ 时，$n!$ 中必然包含因子 $p$，所以 $n! \equiv 0 \pmod p$。

这看起来好像问题解决了？但实际比赛中，几乎不会直接问你 $n! \pmod p$ 的值，因为答案通常是 0。更常见的问题是：

1. 计算组合数 $\binom{n}{k} \pmod p$。
2. 计算 $n! / p^{v_p(n!)} \pmod p$，即从 $n!$ 中除去所有因子 $p$ 后再取模。

这两个问题都需要我们更深入地理解 $n!$ 的结构。

3.1 勒让德定理 (Legendre's Formula)

首先，我们需要知道 $n!$ 中到底包含了多少个因子 $p$。勒让德定理给出了答案。

勒让德定理: 素数 $p$ 在 $n!$ 的素因子分解中出现的次数（即 $p$ 的幂次），记作 $v_p(n!)$，等于：

$$v_p(n!) = \sum_{i=1}^{\infty} \lfloor \frac{n}{p^i} \rfloor = \lfloor \frac{n}{p} \rfloor + \lfloor \frac{n}{p^2} \rfloor + \lfloor \frac{n}{p^3} \rfloor + \dots $$

这个和是有限的，因为当 $p^i > n$ 时，$\lfloor n/p^i \rfloor = 0$。

理解： $\lfloor n/p \rfloor$ 计算了 $1, \dots, n$ 中有多少个 $p$ 的倍数。 $\lfloor n/p^2 \rfloor$ 计算了有多少个 $p^2$ 的倍数。这些数至少贡献了两个因子 $p$，其中一个已经在 $\lfloor n/p \rfloor$ 中计算过，这里补充计算了额外的那个因子 $p$。 以此类推，$\lfloor n/p^k \rfloor$ 计算了 $p^k$ 的倍数贡献的第 $k$ 个因子 $p$。

实现：

#include <bits/stdc++.h>

using namespace std;
using ll = long long;

// 计算 v_p(n!)
ll legendre(ll n, ll p) {
    // 输入: n, p (p为素数)
    // 输出: n! 中因子 p 的个数
    if (n < p) return 0;
    ll cnt = 0;
    while (n > 0) {
        n /= p;
        cnt += n;
    }
    return cnt;
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    
    ll n, p;
    // 题意: 计算 n! 中质因子 p 的数量
    cin >> n >> p; 
    cout << legendre(n, p) << endl;
    
    return 0;
}
// 时间复杂度: O(log_p n)
// 空间复杂度: O(1)

3.2 计算 $n! / p^{v_p(n!)} \pmod p$

这是处理 $n \ge p$ 时最关键的计算之一。我们想计算 $n!$ 中所有不含因子 $p$ 的部分的乘积模 $p$。

令 $f(n) = n! / p^{v_p(n!)}$。我们要求 $f(n) \pmod p$。

考虑 $n! = 1 \times 2 \times \dots \times n$。我们可以把这些数分成两类：

1. 能被 $p$ 整除的数: $p, 2p, 3p, \dots, \lfloor n/p \rfloor p$。
2. 不能被 $p$ 整除的数。

$n! = (\text{不被 } p \text{ 整除的数的乘积}) \times (\text{被 } p \text{ 整除的数的乘积})$

$(\text{被 } p \text{ 整除的数的乘积}) = (p) \times (2p) \times \dots \times (\lfloor n/p \rfloor p)$ $= p^{\lfloor n/p \rfloor} \times (1 \times 2 \times \dots \times \lfloor n/p \rfloor)$ $= p^{\lfloor n/p \rfloor} \times (\lfloor n/p \rfloor)!$

代入 $n!$ 的表达式： $n! = (\text{不被 } p \text{ 整除的数的乘积}) \times p^{\lfloor n/p \rfloor} \times (\lfloor n/p \rfloor)!$

两边同时除以 $p^{v_p(n!)}$： $f(n) = \frac{n!}{p^{v_p(n!)}} = (\text{不被 } p \text{ 整除的数的乘积}) \times \frac{p^{\lfloor n/p \rfloor} \times (\lfloor n/p \rfloor)!}{p^{v_p(n!)}}$

根据勒让德公式 $v_p(n!) = \lfloor n/p \rfloor + v_p(\lfloor n/p \rfloor!)$，我们有 $v_p(n!) - \lfloor n/p \rfloor = v_p(\lfloor n/p \rfloor!)$。 所以 $\frac{p^{\lfloor n/p \rfloor}}{p^{v_p(n!)}} = \frac{1}{p^{v_p(\lfloor n/p \rfloor!)}}$。

代入得到： $f(n) = (\text{不被 } p \text{ 整除的数的乘积}) \times \frac{(\lfloor n/p \rfloor)!}{p^{v_p(\lfloor n/p \rfloor!)}}$ $f(n) = (\text{不被 } p \text{ 整除的数的乘积}) \times f(\lfloor n/p \rfloor)$

现在我们需要计算 $(\text{不被 } p \text{ 整除的数的乘积}) \pmod p$。 考虑 $1, 2, \dots, n$ 中的这些数。我们将它们按模 $p$ 分组。 例如，考虑 $1, \dots, p-1$。它们的乘积是 $(p-1)!$。 考虑 $p+1, \dots, 2p-1$。模 $p$ 后，它们对应 $1, \dots, p-1$。所以它们的乘积模 $p$ 也是 $(p-1)!$。 ... 考虑 $(\lfloor n/p \rfloor - 1)p + 1, \dots, \lfloor n/p \rfloor p - 1$。模 $p$ 后也是 $1, \dots, p-1$。乘积模 $p$ 也是 $(p-1)!$。 最后，还有一部分数：$\lfloor n/p \rfloor p + 1, \dots, n$。模 $p$ 后对应 $1, 2, \dots, n \pmod p$。它们的乘积是 $(n \pmod p)!$。

所以，$(\text{不被 } p \text{ 整除的数的乘积}) \pmod p$ $= [(p-1)!]^{\lfloor n/p \rfloor} \times (n \pmod p)! \pmod p$

根据威尔逊定理，$(p-1)! \equiv -1 \pmod p$。 因此，$(\text{不被 } p \text{ 整除的数的乘积}) \pmod p \equiv (-1)^{\lfloor n/p \rfloor} \times (n \pmod p)! \pmod p$。

将这个结果代回 $f(n)$ 的递推式： $f(n) \pmod p \equiv [(-1)^{\lfloor n/p \rfloor} \times (n \pmod p)!] \times f(\lfloor n/p \rfloor) \pmod p$

这是一个递归计算 $f(n) \pmod p$ 的方法。边界条件是 $f(0) = 1$。

注意： $(-1)^{\lfloor n/p \rfloor}$ 可以写成 $p-1$ 如果 $\lfloor n/p \rfloor$ 是奇数，写成 $1$ 如果 $\lfloor n/p \rfloor$ 是偶数。但有个坑：如果 $p=2$，那么 $-1 \equiv 1 \pmod 2$。所以在 $p=2$ 的情况下，这一项总是 $1$。不过更简单的处理是， $(-1)^{\lfloor n/p \rfloor}$ 在模 $p$ 意义下，只有当 $p > 2$ 且 $\lfloor n/p \rfloor$ 为奇数时才等于 $p-1$，否则等于 $1$。

实现：

#include <bits/stdc++.h>

using namespace std;
using ll = long long;

// 预计算 1! mod p, 2! mod p, ..., (p-1)! mod p
vector<ll> fact;
ll P; // 模数 p

// 快速幂
ll power(ll a, ll b) {
    ll res = 1;
    a %= P;
    while (b > 0) {
        if (b & 1) res = (res * a) % P;
        a = (a * a) % P;
        b >>= 1;
    }
    return res;
}

// 递归计算 n! / p^(v_p(n!)) mod p
ll fact_mod_p_removed(ll n) {
    // 输入: n
    // 输出: n! 中除去所有因子 p 后，模 p 的值
    if (n == 0) return 1;

    // 计算 (n mod p)! mod p
    ll term1 = fact[n % P];

    // 计算 (-1)^(n/p) mod p, 注意 p=2 的情况
    ll term2 = 1;
    if ( (n / P) % 2 == 1 && P > 2) { // 只有 p > 2 且 n/p 是奇数时，才是 -1 (即 p-1)
        term2 = P - 1;
    }
    // 另一种写法避免判断p=2: 
    // if ((n / P) & 1) term2 = P - 1; else term2 = 1;
    // if (P == 2) term2 = 1; // 特判 p=2

    // 递归计算 f(n/p) mod p
    ll term3 = fact_mod_p_removed(n / P);

    return (((term1 * term2) % P) * term3) % P;
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);

    ll n;
    // 题意: 计算 n! / p^(v_p(n!)) mod p，其中 p 是素数
    cin >> n >> P; 

    // 预计算阶乘模 p
    fact.resize(P);
    fact[0] = 1;
    for (int i = 1; i < P; ++i) {
        fact[i] = (fact[i - 1] * i) % P;
    }

    cout << fact_mod_p_removed(n) << endl;

    return 0;
}
// 时间复杂度: O(p + log_p n) (预处理 O(p), 递归 O(log_p n))
// 空间复杂度: O(p) (存储预处理的阶乘)

这个递归方法非常高效，时间复杂度是对数级的。

3.3 卢卡斯定理 (Lucas's Theorem)

卢卡斯定理是计算组合数 $\binom{n}{k} \pmod p$ （$p$ 为素数）的强大武器，尤其在 $n, k$ 很大但 $p$ 相对较小时。

卢卡斯定理: 对于非负整数 $n, k$ 和素数 $p$，有：

$$\binom{n}{k} \equiv \binom{\lfloor n/p \rfloor}{\lfloor k/p \rfloor} \times \binom{n \pmod p}{k \pmod p} \pmod p $$

理解： 这个定理将 $\binom{n}{k} \pmod p$ 的计算分解为两个更小的组合数计算：一个是规模缩小 $p$ 倍的 $\binom{\lfloor n/p \rfloor}{\lfloor k/p \rfloor} \pmod p$，另一个是规模在 $p$ 以内的 $\binom{n \pmod p}{k \pmod p} \pmod p$。

我们可以递归地应用这个定理，直到 $n, k$ 都变成 0。

计算 $\binom{a}{b} \pmod p$ 其中 $a, b < p$： 这需要计算 $a! \times (b!)^{-1} \times ((a-b)! )^{-1} \pmod p$。 由于 $a, b < p$，所以 $b!$ 和 $(a-b)!$ 都不含因子 $p$，它们的模 $p$ 逆元存在。我们可以预计算 $0!, 1!, \dots, (p-1)! \pmod p$，然后用费马小定理（或 exgcd）求逆元。

实现：

#include <bits/stdc++.h>

using namespace std;
using ll = long long;

ll P; // 模数 p

// 快速幂 a^b mod P
ll power(ll a, ll b) {
    ll res = 1;
    a %= P;
    while (b > 0) {
        if (b & 1) res = (res * a) % P;
        a = (a * a) % P;
        b >>= 1;
    }
    return res;
}

// 模逆元 a^-1 mod P (要求 P 是素数)
ll inv(ll a) {
    return power(a, P - 2);
}

// 预计算阶乘和阶乘逆元模 P
vector<ll> fact, invFact;

// 计算 C(n, k) mod P, 其中 n, k < P
ll C_small(ll n, ll k) {
    if (k < 0 || k > n) return 0;
    if (k == 0 || k == n) return 1;
    if (invFact[k] == 0 || invFact[n - k] == 0) { // 逆元不存在（理论上不会发生，因为 k, n-k < P）
         return -1; // Error signal
    }
    // C(n, k) = n! / (k! * (n-k)!) = fact[n] * invFact[k] * invFact[n-k] mod P
    return (((fact[n] * invFact[k]) % P) * invFact[n - k]) % P;
}

// 卢卡斯定理计算 C(n, k) mod P
ll lucas(ll n, ll k) {
    // 输入: n, k, 全局 P (素数)
    // 输出: C(n, k) mod P
    if (k < 0 || k > n) return 0;
    if (k == 0) return 1;
    
    // 递归计算
    ll C1 = C_small(n % P, k % P);
    ll C2 = lucas(n / P, k / P);
    
    return (C1 * C2) % P;
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);

    ll n, k;
    // 题意: 计算 C(n, k) mod p，其中 p 是素数
    cin >> n >> k >> P;

    // 预计算 0! 到 (P-1)! 模 P 及其逆元
    fact.resize(P);
    invFact.resize(P);
    fact[0] = 1;
    invFact[0] = 1; // 0! 的逆元认为是 1
    for (int i = 1; i < P; ++i) {
        fact[i] = (fact[i - 1] * i) % P;
        // 注意: 必须在计算完所有 fact 后再计算 invFact，或者用 exgcd 单独求逆元
        // 这里我们先算完 fact
    }
    // 计算 (P-1)! 的逆元，然后递推计算其他逆元
    invFact[P - 1] = inv(fact[P - 1]);
    for (int i = P - 2; i >= 1; --i) {
        // inv((i)!) = inv((i+1)!) * (i+1) mod P
        invFact[i] = (invFact[i + 1] * (i + 1)) % P;
    }


    cout << lucas(n, k) << endl;

    return 0;
}
// 时间复杂度: 预处理 O(P log P) (费马小定理求逆元) 或 O(P) (递推求逆元), 查询 O(log_p n)
// 空间复杂度: O(P)

卢卡斯定理提供了一种在模素数 $p$ 下计算大组合数的标准方法。

四、 进阶挑战：p 为合数

当模数 $p$ 是合数时，情况变得复杂得多。费马小定理失效，威尔逊定理失效，卢卡斯定理也直接失效。

核心思想是：

1. 将 $p$ 进行素因子分解：$p = p_1^{a_1} p_2^{a_2} \dots p_m^{a_m}$。
2. 分别计算 $n! \pmod{p_i^{a_i}}$ 的值，得到 $m$ 个结果 $r_1, r_2, \dots, r_m$。
3. 利用中国剩余定理 (Chinese Remainder Theorem, CRT) 将这 $m$ 个同余方程：$$x \equiv r_1 \pmod{p_1^{a_1}}$$ $$x \equiv r_2 \pmod{p_2^{a_2}}$$ $$\dots$$ $$x \equiv r_m \pmod{p_m^{a_m}}$$ 合并成一个解 $x \pmod p$。这个 $x$ 就是 $n! \pmod p$。

现在，主要矛盾转化为如何计算 $n! \pmod{p^k}$，其中 $p$ 是素数，$k$ 是正整数。

4.1 计算 $n! \pmod{p^k}$

这可以说是阶乘取模问题中最复杂的部分。我们不能简单地像模素数那样操作，因为模 $p^k$ 意义下的逆元不一定存在（当数是 $p$ 的倍数时）。

同样，我们将 $n!$ 分解： $n! = (\text{不含因子 } p \text{ 的项的乘积}) \times (\text{含有因子 } p \text{ 的项的乘积})$

令 $P = p^k$。 $(\text{含有因子 } p \text{ 的项的乘积}) = (p) \times (2p) \times \dots \times (\lfloor n/p \rfloor p)$ $= p^{\lfloor n/p \rfloor} \times (\lfloor n/p \rfloor)!$

令 $g(n) = \text{不含因子 } p \text{ 的项的乘积} \pmod{P}$。 则 $n! \equiv g(n) \times p^{\lfloor n/p \rfloor} \times (\lfloor n/p \rfloor)! \pmod{P}$ （注意，这里的同余是在模 $P=p^k$ 下讨论，但 $p^{\lfloor n/p \rfloor}$ 可能非常大，这里只是形式化表示）

我们真正需要的是计算 $n! \pmod{p^k}$ 的值。直接的递推似乎很困难。

一个更有效的方法是分别计算：

1. $v_p(n!)$：$n!$ 中因子 $p$ 的个数（用勒让德公式）。
2. $n! / p^{v_p(n!)} \pmod{p^k}$：$n!$ 中除去所有因子 $p$ 后，模 $p^k$ 的值。

令 $F(n, p, k) = n! / p^{v_p(n!)} \pmod{p^k}$。 我们有递推关系： $n! = (\prod_{i=1, p \nmid i}^{n} i) \times (p \cdot 2p \cdot \dots \cdot \lfloor n/p \rfloor p)$ $n! = (\prod_{i=1, p \nmid i}^{n} i) \times p^{\lfloor n/p \rfloor} \times (\lfloor n/p \rfloor)!$

两边同除以 $p^{v_p(n!)}$，并利用 $v_p(n!) = \lfloor n/p \rfloor + v_p(\lfloor n/p \rfloor!)$： $\frac{n!}{p^{v_p(n!)}} = (\prod_{i=1, p \nmid i}^{n} i) \times \frac{p^{\lfloor n/p \rfloor}}{p^{v_p(n!)}} \times (\lfloor n/p \rfloor)!$ $F(n, p, k) \equiv (\prod_{i=1, p \nmid i}^{n} i) \times \frac{(\lfloor n/p \rfloor)!}{p^{v_p(\lfloor n/p \rfloor!)}} \pmod{p^k}$ $F(n, p, k) \equiv (\prod_{i=1, p \nmid i}^{n} i) \times F(\lfloor n/p \rfloor, p, k) \pmod{p^k}$

关键在于计算 $\prod_{i=1, p \nmid i}^{n} i \pmod{p^k}$。 我们将 $1, \dots, n$ 分成若干个长度为 $P=p^k$ 的块。 考虑一个完整的块 $[mP+1, (m+1)P]$。其中与 $p$ 互素的数的乘积模 $P$ 是多少？ $\prod_{i=mP+1, p \nmid i}^{(m+1)P} i \equiv \prod_{j=1, p \nmid j}^{P} (mP+j) \pmod P$ $\equiv \prod_{j=1, p \nmid j}^{P} j \pmod P$ 这说明，在每个长度为 $P$ 的完整区间内，与 $p$ 互素的数的乘积模 $P$ 是相同的！令这个值为 $Val = \prod_{j=1, p \nmid j}^{P} j \pmod P$。

那么 $\prod_{i=1, p \nmid i}^{n} i$ 可以写成： $= (\prod_{j=1, p \nmid j}^{P} j)^{\lfloor n/P \rfloor} \times (\prod_{i=\lfloor n/P \rfloor P + 1, p \nmid i}^{n} i) \pmod P$ $= Val^{\lfloor n/P \rfloor} \times (\prod_{j=1, p \nmid j}^{n \pmod P} j) \pmod P$

现在，我们需要计算两个部分：

1. $Val = \prod_{j=1, p \nmid j}^{p^k} j \pmod{p^k}$。
2. $Term(m) = \prod_{j=1, p \nmid j}^{m} j \pmod{p^k}$ 对于 $m < p^k$。

计算 $Term(m)$ 可以通过暴力 $O(m)$ 完成。 计算 $Val$ 也可以暴力 $O(p^k)$ 完成。

结合起来，计算 $F(n, p, k)$ 的递归函数：

// 假设 P = p^k
ll calculate_F(ll n, ll p, ll P) {
    if (n == 0) return 1;

    // 计算 Term(n % P) = (product_{j=1, p not dividing j}^{n mod P} j) mod P
    ll term_last_block = 1;
    for (ll i = 1; i <= n % P; ++i) {
        if (i % p != 0) {
            term_last_block = (term_last_block * i) % P;
        }
    }

    // 计算 Val = (product_{j=1, p not dividing j}^{P} j) mod P
    // 这个 Val 可以在递归外预计算一次，或者每次递归都算？
    // 实际上 Val 在模 p^k 下有个性质，可以优化，但我们先用暴力的思想
    ll val_full_block = 1;
    // 优化：由于 (m*p^k + j) mod p^k = j，所以完整的块的乘积在模 p^k 下都一样
    // 且 Val = product_{j=1, p not dividing j}^{p^k} j mod p^k
    // 这个值有个性质：若 p>2 或 k<3, Val = -1 mod p^k; 否则 Val = 1 mod p^k
    // 我们这里先写成需要计算的样子
    if (n >= P) { // 只有当 n >= P 时才需要计算完整块
        for (ll i = 1; i <= P; ++i) {
            if (i % p != 0) {
                val_full_block = (val_full_block * i) % P;
            }
        }
    } else {
         val_full_block = 1; // 如果 n < P，没有完整的块
    }
    
    // 计算 Val^(n/P) mod P
    ll full_blocks_prod = power(val_full_block, n / P, P); // 需要快速幂

    // 递归计算 F(n/p, p, k)
    ll recursive_part = calculate_F(n / p, p, P);

    // 最终结果 F(n) = (Val^(n/P) * Term(n%P)) * F(n/p) mod P
    ll res = (full_blocks_prod * term_last_block) % P;
    res = (res * recursive_part) % P;
    return res;
}

// 然后 n! mod P = (F(n, p, k) * p^(v_p(n!)) ) mod P ?
// 不对，直接这样乘 p^(v_p(n!)) 会出问题，因为 F(n, p, k) 是模 P 计算的
// 而且 p^(v_p(n!)) 可能 >= P

正确的处理 $n! \pmod{p^k}$ 的方法：

1. 计算 $v_p(n!)$ 使用勒让德公式。记为 vp_n.

2. 如果 $v_p(n!) \ge k$，那么 $p^k$ 是 $n!$ 的因子， $n! \equiv 0 \pmod{p^k}$。

3. 如果 $v_p(n!) < k$，我们需要计算 $n! \pmod{p^k}$。这才是难点。 我们需要计算 $F(n, p, k) = n! / p^{v_p(n!)} \pmod{p^k}$。 然后最终结果是 $(F(n, p, k) \times p^{v_p(n!)}) \pmod{p^k}$。这里的乘法可以直接做，因为 $v_p(n!) < k$，所以 $p^{v_p(n!)}$ 不会引入模 $p^k$ 的问题（即不会变成 0）。

   如何计算 $F(n, p, k) \pmod{p^k}$？ 上面推导的递推关系是正确的： $F(n, p, k) \equiv (\prod_{i=1, p \nmid i}^{n} i) \times F(\lfloor n/p \rfloor, p, k) \pmod{p^k}$

   计算 $\prod_{i=1, p \nmid i}^{n} i \pmod{p^k}$： 令 $P = p^k$. $= (\prod_{j=1, p \nmid j}^{P} j)^{\lfloor n/P \rfloor} \times (\prod_{j=1, p \nmid j}^{n \pmod P} j) \pmod P$

   这个方法是可行的，但每次递归都需要 $O(P)$ 的计算（如果 $Val$ 和 $Term(m)$ 都暴力算）。总复杂度会是 $O(P \log_p n)$。如果 $P=p^k$ 很大，这不行。

   优化 $\prod_{i=1, p \nmid i}^{n} i \pmod{p^k}$ 的计算: 我们把 $1, \dots, n$ 写成 $n = qp + r$ ($0 \le r < p$)。 考虑 $1, \dots, n$ 中的数，按模 $p$ 分类。

   • $1, 2, \dots, p-1$
   • $p+1, \dots, 2p-1$
   • ...
   • $(q-1)p+1, \dots, qp-1$
   • $qp+1, \dots, qp+r = n$

   每一组完整的 $(m-1)p+1, \dots, mp-1$ 的乘积模 $p$ 是 $(p-1)! \equiv -1 \pmod p$。 但我们现在需要模 $p^k$。

   一个更标准的做法是（以计算 $\binom{n}{m} \pmod{p^k}$ 为例，这涉及 $n!, m!, (n-m)! \pmod{p^k}$ 的计算）：

   扩展卢卡斯定理 (Extended Lucas Theorem) 的核心步骤

   目标：计算 $\binom{n}{m} \pmod{p^k}$. $\binom{n}{m} = \frac{n!}{m!(n-m)!}$

   1. 计算 $v_p(n!), v_p(m!), v_p((n-m)!)$。 记为 $v_n, v_m, v_{nm}$。
   2. 计算 $\binom{n}{m}$ 中因子 $p$ 的总幂次 $E = v_n - v_m - v_{nm}$。
   3. 计算 $N = n! / p^{v_n} \pmod{p^k}$。
   4. 计算 $M = m! / p^{v_m} \pmod{p^k}$。
   5. 计算 $NM = (n-m)! / p^{v_{nm}} \pmod{p^k}$。
   6. 计算 $M^{-1} \pmod{p^k}$ 和 $NM^{-1} \pmod{p^k}$。这需要 $M, NM$ 与 $p$ 互素，这是由定义保证的。需要用 exgcd 求逆元。
   7. 最终结果是 $\binom{n}{m} \equiv N \times M^{-1} \times NM^{-1} \times p^E \pmod{p^k}$。

   核心就是计算 $F(x, p, k) = x! / p^{v_p(x!)} \pmod{p^k}$。

   计算 $F(n, p, k) \pmod{p^k}$ 的高效方法

   令 $P = p^k$。 $F(n) \equiv (\prod_{i=1, p \nmid i}^{n} i) \times F(\lfloor n/p \rfloor) \pmod P$

   $\prod_{i=1, p \nmid i}^{n} i = (\prod_{j=1, p \nmid j}^{P} j)^{\lfloor n/P \rfloor} \times (\prod_{j=1, p \nmid j}^{n \pmod P} j) \pmod P$

   • $\prod_{j=1, p \nmid j}^{n \pmod P} j \pmod P$：可以 $O(P)$ 预处理 $Term(m) = \prod_{j=1, p \nmid j}^{m} j \pmod P$ for $m=1 \dots P-1$。查询 O(1)。
   • $\prod_{j=1, p \nmid j}^{P} j \pmod P$：可以 $O(P)$ 预处理一次得到 $Val$。

   递归计算 $F(n)$ 的复杂度将是 $O(P + \log_p n)$。

   实现扩展卢卡斯需要的函数：

   1. power(a, b, mod): 快速幂
   2. exgcd(a, b, &x, &y): 扩展欧几里得算法
   3. inv(a, mod): 使用 exgcd 求模 mod 的逆元
   4. legendre(n, p): 计算 $v_p(n!)$
   5. calc_F(n, p, pk): 计算 $n! / p^{v_p(n!)} \pmod{p^k}$
   6. C_mod_pk(n, m, p, pk): 计算 $\binom{n}{m} \pmod{p^k}$
   7. CRT(remainders, moduli): 中国剩余定理合并结果
   8. ExtendedLucas(n, m, P): 主函数，分解 $P$，调用 C_mod_pk 和 CRT。

   calc_F(n, p, pk) 的实现:

   #include <bits/stdc++.h>
   
   using namespace std;
   using ll = long long;
   
   // 注意：这里的 pk = p^k (模数)
   
   // 快速幂 (模数是 pk)
   ll power(ll a, ll b, ll pk) {
       ll res = 1;
       a %= pk;
       while (b > 0) {
           if (b & 1) res = (res * a) % pk;
           a = (a * a) % pk;
           b >>= 1;
       }
       return res;
   }
   
   // 预计算 Term(m) = product_{j=1, p not dividing j}^{m} j mod pk
   // 预计算 Val = product_{j=1, p not dividing j}^{pk} j mod pk
   vector<ll> term_prod;
   ll val_full_block;
   ll P_prime; // 质数 p
   
   // 计算 n! / p^(v_p(n!)) mod pk
   ll calc_F(ll n, ll pk) {
       if (n == 0) return 1;
   
       // 计算 (product_{j=1, p not dividing j}^{n mod pk} j) mod pk
       ll last_block_prod = term_prod[n % pk];
   
       // 计算 (Val)^(n/pk) mod pk
       ll full_blocks_prod = power(val_full_block, n / pk, pk);
   
       // 递归 F(n/p)
       ll recursive_part = calc_F(n / P_prime, pk); // 注意这里是 n / p
   
       // F(n) = (full_blocks_prod * last_block_prod) * F(n/p) mod pk
       ll res = (full_blocks_prod * last_block_prod) % pk;
       res = (res * recursive_part) % pk;
       return res;
   }
   
   // 预处理 Term 和 Val
   void precompute(ll p, ll pk) {
       P_prime = p;
       term_prod.resize(pk + 1);
       term_prod[0] = 1;
       for (ll i = 1; i <= pk; ++i) {
           if (i % p != 0) {
               term_prod[i] = (term_prod[i - 1] * i) % pk;
           } else {
               term_prod[i] = term_prod[i - 1];
           }
       }
       val_full_block = term_prod[pk];
   }
   

// --- 下面是组合数和 CRT 的部分，为了完整性也给出框架 ---

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

// Modular Inverse using Exgcd
ll mod_inv(ll a, ll mod) {
    ll x, y;
    ll d = exgcd(a, mod, x, y);
    if (d == 1) {
        return (x % mod + mod) % mod;
    } else {
        return -1; // Inverse does not exist
    }
}

// Legendre's Formula
ll legendre(ll n, ll p) {
    if (n < p) return 0;
    ll cnt = 0;
    while (n > 0) {
        n /= p;
        cnt += n;
    }
    return cnt;
}

// Calculate C(n, m) mod pk
ll C_mod_pk(ll n, ll m, ll p, ll pk) {
    if (m < 0 || m > n) return 0;
    if (m == 0 || m == n) return 1;

    precompute(p, pk); // 预处理 Term 和 Val

    ll vn = legendre(n, p);
    ll vm = legendre(m, p);
    ll vnm = legendre(n - m, p);

    ll E = vn - vm - vnm; // p 的幂次

    // 如果 p^E >= pk，结果可能为 0，但这里 E 不一定 >= k
    // 我们需要计算 p^E mod pk
    ll p_pow_E = power(p, E, pk); // 计算 p^E mod pk

    ll fn = calc_F(n, pk);
    ll fm = calc_F(m, pk);
    ll fnm = calc_F(n - m, pk);

    ll inv_fm = mod_inv(fm, pk);
    ll inv_fnm = mod_inv(fnm, pk);

    if (inv_fm == -1 || inv_fnm == -1) {
         // 理论上不应发生，因为 F(x) 保证与 p 互质
         return -1; // Error signal
    }

    ll res = fn;
    res = (res * inv_fm) % pk;
    res = (res * inv_fnm) % pk;
    res = (res * p_pow_E) % pk; // 最后乘上 p^E

    return res;
}

// 中国剩余定理 (CRT) - 求解 x = r_i (mod m_i)
// 返回解 x mod M，其中 M = product(m_i)
ll CRT(const vector<ll>& r, const vector<ll>& m) {
    ll M = 1;
    for (ll mod : m) {
        M *= mod;
    }
    ll result = 0;
    for (size_t i = 0; i < r.size(); ++i) {
        ll Mi = M / m[i];
        ll inv_Mi = mod_inv(Mi, m[i]); // Mi 在模 m[i] 下的逆元
        if (inv_Mi == -1) return -1; // Error if inverse doesn't exist
        result = (result + r[i] * Mi % M * inv_Mi % M) % M;
    }
    return (result + M) % M; // Ensure positive result
}

// 扩展卢卡斯主函数
ll extended_lucas(ll n, ll m, ll P_mod) {
     if (m < 0 || m > n) return 0;
     if (m == 0 || m == n) return 1 % P_mod; // 注意取模
     if (P_mod == 1) return 0;

     vector<pair<ll, ll>> factors; // 存储 P_mod 的质因子分解 {p, p^k}
     ll temp_P = P_mod;
     for (ll i = 2; i * i <= temp_P; ++i) {
         if (temp_P % i == 0) {
             ll pk = 1;
             ll p = i;
             while (temp_P % i == 0) {
                 pk *= i;
                 temp_P /= i;
             }
             factors.push_back({p, pk});
         }
     }
     if (temp_P > 1) { // 剩余一个大的质因子
         factors.push_back({temp_P, temp_P});
     }

     vector<ll> remainders;
     vector<ll> moduli;

     for (const auto& factor : factors) {
         ll p = factor.first;
         ll pk = factor.second;
         ll rem = C_mod_pk(n, m, p, pk);
         if (rem == -1) return -1; // Error in calculation
         remainders.push_back(rem);
         moduli.push_back(pk);
     }
     
     // 如果只有一个因子，CRT 不需要
     if (remainders.size() == 1) {
         return remainders[0];
     }

     return CRT(remainders, moduli);
}

int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);

    ll n, m, P_mod;
    // 题意: 计算 C(n, m) mod P_mod, P_mod 可以是合数
    cin >> n >> m >> P_mod;

    cout << extended_lucas(n, m, P_mod) << endl;

    return 0;
}
// 时间复杂度: 
// 质因数分解 O(sqrt(P_mod))
// 对每个质因子 p^k:
//   预处理 O(p^k)
//   C_mod_pk 计算 O(p^k + log_p n * log pk) (快速幂和递归调用)
//   求逆元 O(log pk) (exgcd)
// CRT O(num_factors * log P_mod)
// 主要瓶颈在于预处理的 O(p^k)。如果 P_mod 的某个质因子 p 对应的 p^k 很大，会超时。
// 空间复杂度: O(max(p^k)) (预处理 term_prod)
```

4.2 小结与注意事项 (合数模)

• 核心: 分解模数 $P$，计算模 $p^k$ 下的组合数，用 CRT 合并。
• 难点: 计算 $n! / p^{v_p(n!)} \pmod{p^k}$。
• 复杂度: 扩展卢卡斯的主要瓶颈在于计算模 $p^k$ 的部分，特别是预处理阶乘相关项，复杂度与 $p^k$ 相关。如果 $p^k$ 过大（例如 $> 10^6$ 或 $10^7$），则该算法可能超时。
• 逆元: 模 $p^k$ 求逆元必须用 exgcd。
• 精度: 全程使用 long long 防止溢出。

五、 实战技巧与坑点

1. 看清模数性质: 第一时间判断模数 $p$ 是素数还是合数。决定了你使用 Lucas 还是 Extended Lucas。
2. 预处理: 对于素数模 $p$，预处理 $1!, \dots, (p-1)!$ 及其逆元是标准操作，复杂度 O(p)。对于扩展 Lucas，预处理 $\prod_{j=1, p \nmid j}^m j \pmod{p^k}$ 是关键，复杂度 O(p^k)。
3. 边界情况: $n=0, m=0, m>n$ 等边界要正确处理。 $0! = 1$。$\binom{n}{m}=0$ if $m<0$ or $m>n$.
4. long long: 几乎所有中间计算都要用 long long。
5. 模运算时机: 每次乘法后都要及时取模，防止溢出。加法后也要 (a + b) % mod，减法后要 (a - b % mod + mod) % mod 防止负数。
6. 复杂度意识: 估算 $n, k, p$ (或 $p^k$) 的大小，判断 O(p), O(log n), O(p^k) 等复杂度是否可以接受。
7. 威尔逊定理的应用场景: 虽然不直接计算阶乘，但它在推导 $n! / p^{v_p(n!)} \pmod p$ 时用到 $(p-1)! \equiv -1 \pmod p$。某些特殊题目可能直接考察威尔逊定理。

六、 总结

我们系统地梳理了阶乘取模的几种情况：

1. $n < p$, $p$ 是素数: 直接计算 $O(n)$。
2. $n \ge p$, $p$ 是素数:
   • $n! \equiv 0 \pmod p$。
   • 计算 $n! / p^{v_p(n!)} \pmod p$：使用基于威尔逊定理的递归方法 $O(p + \log_p n)$。
   • 计算 $\binom{n}{k} \pmod p$：使用卢卡斯定理 $O(p + \log_p n)$。
3. $p$ 是合数:
   • 分解 $p = p_1^{a_1} \dots p_m^{a_m}$。
   • 计算 $\binom{n}{k} \pmod{p_i^{a_i}}$ (扩展卢卡斯的核心步骤)：
     • 需要计算 $x! / p^{v_p(x!)} \pmod{p^k}$。
     • 复杂度依赖 $p^k$ 的大小，预处理 $O(p^k)$。
   • 使用 CRT 合并结果。