洛谷 P3807 【模板】卢卡斯定理 / Lucas 定理

// ✗ 第 ② 版：把[本章第 6 步](/ch/43-combinatorics/)那套模板原样搬过来 —— 而它在这道题上**恒输出 0**
//
// ★★★ 这道题就是本章第 6 步那句警告的**考场现场**：
//   那套「阶乘 + 逆元」要求 **p > n**。而这里 p ≤ 10⁵、n+m 可以到 2×10⁵
//   ⇒ 只要 **p ≤ n+m**，(n+m)! 里就含着因子 p ⇒ fac[n+m] ≡ 0 (mod p)
//   ⇒ 答案 = 0 · 逆元 · 逆元 = **0**，而且是**恒等于 0**，不是「有时错」。
//
// ⚠ 更糟的是逆元那一步：fac[n] 也可能 ≡ 0，而 **0 在模 p 下没有逆元** ——
//   qpow(0, p−2) 老老实实算出来是 0，程序一声不吭，答案却早就没有意义了。
//
// ⇒ 这一版**能跑、不崩、输出格式完全正常**，它只是**一直在答 0**。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

static ll P;
static ll fac[200005], inv[200005];

static ll qpow(ll a, ll b) {
    ll res = 1 % P; a %= P;
    while (b) { if (b & 1) res = res * a % P; a = a * a % P; b >>= 1; }
    return res;
}

int main() {
    int T;
    if (scanf("%d", &T) != 1) return 0;
    while (T--) {
        ll n, m;
        if (scanf("%lld %lld %lld", &n, &m, &P) != 3) return 0;
        ll N = n + m;
        fac[0] = 1 % P;
        for (ll i = 1; i <= N; i++) fac[i] = fac[i - 1] * i % P;
        inv[N] = qpow(fac[N], P - 2);
        for (ll i = N; i >= 1; i--) inv[i - 1] = inv[i] * i % P;
        printf("%lld\n", fac[N] * inv[n] % P * inv[m] % P);
    }
    return 0;
}

// P3807 【模板】卢卡斯定理 —— 正解
//
// 要算 C(n+m, n) mod p，题面保证 p 是质数，而 **p ≤ 10⁵、n+m 可以到 2×10⁵**
// ⇒ [本章第 6 步](/ch/43-combinatorics/)那套「阶乘 + 逆元」的前提（**p > n**）**在这里不成立**：
//   (n+m)! 里含着因子 p ⇒ 它 ≡ 0 (mod p) ⇒ 那一版会恒输出 0（见 p3807Fac.cpp）。
//
// ★ Lucas 定理：把 n 和 m 按 **p 进制**拆开，逐位求组合数再乘起来：
//     C(n, m) ≡ C(n mod p, m mod p) · C(n / p, m / p)   (mod p)
//   每一位都 < p ⇒ **每一位都退回到了「p > n」那个成立的情形**，可以用阶乘 + 逆元。
//
// ⚠ 三处容易漏：
//   ① 每组数据的 p 不一样 ⇒ **阶乘表要重算**（见 p3807Once.cpp）；
//   ② 递归那一半不能少（见 p3807OneDigit.cpp）；
//   ③ 单位那一层里 m 可能比 n 大 ⇒ **必须返回 0**（见 p3807NoGuard.cpp）。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

static ll P;
static ll fac[100005], inv[100005];

static ll qpow(ll a, ll b) {
    ll res = 1 % P; a %= P;
    while (b) { if (b & 1) res = res * a % P; a = a * a % P; b >>= 1; }
    return res;
}

/** 只到 p−1 就够了：Lucas 每一层的两个下标都 < p */
static void prep() {
    fac[0] = 1 % P;
    for (ll i = 1; i < P; i++) fac[i] = fac[i - 1] * i % P;
    inv[P - 1] = qpow(fac[P - 1], P - 2);
    for (ll i = P - 1; i >= 1; i--) inv[i - 1] = inv[i] * i % P;
}

/** 单位组合数：这里 n、m 都 < p，前提成立 */
static ll C(ll n, ll m) {
    if (m > n) return 0;                       // ★ ③ 这一行少不得
    return fac[n] * inv[m] % P * inv[n - m] % P;
}

static ll lucas(ll n, ll m) {
    if (m == 0) return 1 % P;
    return C(n % P, m % P) * lucas(n / P, m / P) % P;   // ★ ② 递归那一半
}

int main() {
    int T;
    if (scanf("%d", &T) != 1) return 0;
    while (T--) {
        ll n, m;
        if (scanf("%lld %lld %lld", &n, &m, &P) != 3) return 0;
        prep();                                // ★ ① 每组都要重算
        printf("%lld\n", lucas(n + m, n));
    }
    return 0;
}

// P3807 的度量程序：./p3807Count csv        （本页的数字都出自它）
//
//   kummer: ★★★ 一条**和 Lucas 一行代码都不共享**的路 —— Kummer 定理：
//           C(n+m, n) ≡ 0 (mod p)  ⟺  n 和 m 按 **p 进制**相加时**有进位**。
//           拿它把正解验一遍（小范围全枚举），再用它说清楚 NoGuard 那个错法的抓获率。
//   trig  : ★ 四个错法各自的**触发条件**在各档里出现几轮。
//   cost  : ★ 三种写法的规模账（杨辉三角的格子数、两种阶乘表的长度、顶格秒表）。
//
// ⚠ 这里**复刻**了 p3807Gen.cpp 的档位逻辑（同一个 mt19937、同一个种子公式）。
// ⚠ 秒表用 steady_clock 在进程内量，每档 3 次取中位数
//   （[第 12 章 P1923 那一跤](/sol/p1923/)：秒表断言的主语是余量，不是量级）。
#include <bits/stdc++.h>
#include <chrono>
using namespace std;
using namespace std::chrono;
typedef long long ll;

static double med3(double a, double b, double c) { return max(min(a, b), min(max(a, b), c)); }

static vector<int> primesUpTo(int n) {
    vector<char> isc(n + 1, 0);
    vector<int> ps;
    for (int i = 2; i <= n; i++) {
        if (!isc[i]) ps.push_back(i);
        for (ll j = (ll)i * i; j <= n; j += i) isc[j] = 1;
    }
    return ps;
}

/** 复刻 p3807Gen.cpp：把一轮的 T 组数据吐出来 */
static void gen(unsigned seed, int mode, vector<array<ll, 3>>& out) {
    mt19937 rng(seed * 1000003u + 20260903u);
    int nmMax = mode == 0 ? 20 : (mode == 3 ? 100000 : 1500);
    int pMax  = mode == 0 ? 100 : (mode == 3 ? 100000 : 200);
    vector<int> ps = primesUpTo(pMax);
    int T = 1 + (int)(rng() % 10u);
    unsigned pick = rng() % (unsigned)ps.size();
    int shared = ps[pick];
    out.clear();
    for (int t = 0; t < T; t++) {
        ll n, m, p;
        if (mode == 2) {
            unsigned lo = (unsigned)(ps.size() / 2);
            unsigned q = lo + rng() % (unsigned)(ps.size() - lo);
            p = ps[q];
            unsigned r = rng() % (unsigned)p;
            ll N = (ll)r;
            unsigned r2 = rng() % (unsigned)(N + 1);
            n = (ll)r2; m = N - n;
        } else {
            unsigned r1 = rng() % (unsigned)nmMax, r2 = rng() % (unsigned)nmMax;
            n = 1 + (ll)r1; m = 1 + (ll)r2;
            if (mode == 0) p = shared;
            else { unsigned r3 = rng() % (unsigned)ps.size(); p = ps[r3]; }
        }
        out.push_back({n, m, p});
    }
}

/** Kummer：n 和 m 按 p 进制相加，有没有进位 */
static bool hasCarry(ll n, ll m, ll p) {
    ll carry = 0;
    while (n || m) {
        ll d = n % p + m % p + carry;
        carry = d >= p ? 1 : 0;
        if (carry) return true;
        n /= p; m /= p;
    }
    return false;
}

static ll P;
static vector<ll> fac, inv;
static ll qpow(ll a, ll b) {
    ll res = 1 % P; a %= P;
    while (b) { if (b & 1) res = res * a % P; a = a * a % P; b >>= 1; }
    return res;
}
static void prep() {
    fac.assign((size_t)P, 0); inv.assign((size_t)P, 0);
    fac[0] = 1 % P;
    for (ll i = 1; i < P; i++) fac[i] = fac[i - 1] * i % P;
    inv[P - 1] = qpow(fac[P - 1], P - 2);
    for (ll i = P - 1; i >= 1; i--) inv[i - 1] = inv[i] * i % P;
}
static ll Cc(ll n, ll m) {
    if (m > n) return 0;
    return fac[n] * inv[m] % P * inv[n - m] % P;
}
static ll lucas(ll n, ll m) {
    if (m == 0) return 1 % P;
    return Cc(n % P, m % P) * lucas(n / P, m / P) % P;
}

/** 杨辉三角（滚动两行），只在小 N 上用 */
static ll pascalC(ll N, ll k, ll p) {
    vector<ll> prev(N + 1, 0), cur(N + 1, 0);
    prev[0] = 1 % p;
    for (ll i = 1; i <= N; i++) {
        cur[0] = 1 % p;
        for (ll j = 1; j <= i; j++) cur[j] = (prev[j - 1] + prev[j]) % p;
        swap(prev, cur);
    }
    return prev[k];
}

int main(int argc, char** argv) {
    bool csv = (argc > 1 && string(argv[1]) == "csv");

    /* ① Kummer 校验：小范围全枚举，「答案 ≡ 0」和「有进位」对不对得上 */
    int kumCases = 0, kumZero = 0, kumCarry = 0, kumBad = 0, pascalBad = 0;
    for (int p : {2, 3, 5, 7, 11, 13})
        for (ll n = 0; n <= 40; n++)
            for (ll m = 0; m <= 40; m++) {
                P = p; prep();
                ll a = lucas(n + m, n);
                ll b = pascalC(n + m, n, p);
                bool carry = hasCarry(n, m, (ll)p);
                kumCases++;
                if (a == 0) kumZero++;
                if (carry) kumCarry++;
                if ((a == 0) != carry) kumBad++;
                if (a != b) pascalBad++;
            }

    /* ② 四个档 × 四个触发条件 */
    int trigFac[4] = {0, 0, 0, 0}, trigOnce[4] = {0, 0, 0, 0}, trigCarry[4] = {0, 0, 0, 0};
    int facBad[4] = {0, 0, 0, 0}, oneBad[4] = {0, 0, 0, 0};
    for (int mode = 0; mode < 4; mode++) {
        vector<array<ll, 3>> rows;
        for (int s = 1; s <= 300; s++) {
            gen((unsigned)s, mode, rows);
            bool fac_ = false, once_ = false, carry_ = false, facB = false, oneB = false;
            for (auto& r : rows) if (r[2] <= r[0] + r[1]) fac_ = true;
            for (size_t i = 1; i < rows.size(); i++) if (rows[i][2] != rows[0][2]) once_ = true;
            for (auto& r : rows) if (hasCarry(r[0], r[1], r[2])) carry_ = true;
            for (auto& r : rows) {
                ll n = r[0], m = r[1];
                P = r[2]; prep();
                ll N = n + m;
                ll truth = lucas(N, n);
                // 模板版：p ≤ N 时 fac[N] ≡ 0 ⇒ 它恒输出 0
                ll facAns = (P <= N) ? 0 : truth;
                if (facAns != truth) facB = true;
                // 只算一位那版：漏掉了 lucas(N / p, n / p) 那一乘
                ll oneAns = Cc(N % P, n % P);
                if (oneAns != truth) oneB = true;
            }
            if (fac_) trigFac[mode]++;
            if (once_) trigOnce[mode]++;
            if (carry_) trigCarry[mode]++;
            if (facB) facBad[mode]++;
            if (oneB) oneBad[mode]++;
        }
    }

    /* ③ 规模账 */
    const ll NM = 200000;                       // n + m 顶格
    double gbPascalOps = (double)NM * (double)NM;
    ll facLenTemplate = NM + 1;                 // 模板版的阶乘表要开到 n+m
    ll facLenLucas = 100000;                    // Lucas 只要开到 p

    double msLucas;
    {
        double t[3];
        vector<int> ps = primesUpTo(100000);
        for (int r = 0; r < 3; r++) {
            auto t0 = steady_clock::now();
            volatile ll sink = 0;
            for (int g = 0; g < 10; g++) {      // 题面顶格：10 组
                P = ps[ps.size() - 1 - g];      // 接近 10⁵ 的质数
                prep();
                sink += lucas(200000, 100000);
            }
            t[r] = duration<double, milli>(steady_clock::now() - t0).count();
        }
        msLucas = med3(t[0], t[1], t[2]);
    }

    if (csv) {
        printf("kumCases,%d\n", kumCases);
        printf("kumZero,%d\n", kumZero);
        printf("kumCarry,%d\n", kumCarry);
        printf("kumBad,%d\n", kumBad);
        printf("pascalBad,%d\n", pascalBad);
        for (int i = 0; i < 4; i++) printf("trigFac%d,%d\n", i, trigFac[i]);
        for (int i = 0; i < 4; i++) printf("trigOnce%d,%d\n", i, trigOnce[i]);
        for (int i = 0; i < 4; i++) printf("trigCarry%d,%d\n", i, trigCarry[i]);
        for (int i = 0; i < 4; i++) printf("facBad%d,%d\n", i, facBad[i]);
        for (int i = 0; i < 4; i++) printf("oneBad%d,%d\n", i, oneBad[i]);
        printf("pascalOps,%.3g\n", gbPascalOps);
        printf("facLenTemplate,%lld\n", facLenTemplate);
        printf("facLenLucas,%lld\n", facLenLucas);
        printf("msLucas,%.1f\n", msLucas);
        printf("lucasUnderLimit,%.0f\n", 1000.0 / msLucas);
        return 0;
    }

    printf("① Kummer 定理：C(n+m, n) ≡ 0 (mod p)  ⟺  n 和 m 按 p 进制相加有进位\n");
    printf("   六个质数 × n, m ∈ [0, 40] 全枚举，共 %d 组：答案为 0 的 %d 组、有进位的 %d 组\n",
           kumCases, kumZero, kumCarry);
    printf("   ⇒ 两边对不上的：**%d 组**（Lucas 和杨辉三角对不上的：%d 组）\n", kumBad, pascalBad);
    printf("   ⇒ ★★ 这条路和 Lucas 一行代码都不共享 —— 它同时验了正解，也解释了 NoGuard 的抓获率\n");

    const char* NAME[4] = {"顺手写的（n,m ≤ 20、全组共用一个 p）", "主力档（n,m ≤ 1500、每组一个 p）",
                           "⚠ p > n + m", "顶格（n,m ≤ 10⁵）"};
    printf("\n② 四个档各 300 轮，四个触发条件各出现几轮\n");
    for (int i = 0; i < 4; i++)
        printf("   档 %d %-38s p ≤ n+m 的 %3d 轮（其中真会答错 %3d）／ 各组 p 不全同的 %3d 轮 ／ 有进位的 %3d 轮（只算一位真会错 %3d）\n",
               i, NAME[i], trigFac[i], facBad[i], trigOnce[i], trigCarry[i], oneBad[i]);
    printf("   ⇒ ★ 档 0 那个 0（各组 p 全一样）就是「阶乘表只算一次」一次都抓不到的原因\n");
    printf("   ⇒ ★ 档 2 那三个 0 是**一个对照档同时给三个错法做的自检**\n");

    printf("\n③ 规模账（题面顶格：n + m = 2×10⁵、p ≤ 10⁵、T = 10）\n");
    printf("   杨辉三角      要做 %.3g 次加法（一组！题面最多 10 组）  ⇒ 只能当参照物\n", gbPascalOps);
    printf("   本章模板      阶乘表开到 %lld 项 —— 跑得动，可 p ≤ n+m 时**恒输出 0**\n", facLenTemplate);
    printf("   ★ Lucas       阶乘表只开到 %lld 项（p 那么大就够），10 组顶格 %.1f ms ⇒ 余量 %.0f 倍\n",
           facLenLucas, msLucas, 1000.0 / msLucas);
    return 0;
}