洛谷 P1313 [NOIP 2011 提高组] 计算系数

// ✗ 第 ① 版：所有人真实的第一反应 —— 组合数照定义算 C(k,n) = k! / (n! · m!)
//
// 幂那两下大家都会记得取模（题面明写了要模 10007），栽的是**阶乘**这一头：
// long long 装得下 **20!**（2.43×10¹⁸ < 9.22×10¹⁸），**21! 就是第一个装不下的**。
//
// ★ 所以它不是「一交就零分」—— 题面 30% 那一档写着 `0 ≤ k ≤ 10`，
//   那一整档它算得**一个字不差**。⇒ **这一版稳拿 30 分。**
// ⚠ 到 100% 档 k ≤ 1000 时，k! 有 2568 位，long long 早在第 21 步就绕回去了；
//   而绕回去之后那句整数除法连「除得尽」都不成立 —— 输出的是个看着很正常的垃圾数。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll P = 10007;

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

int main() {
    ll a, b, k, n, m;
    if (scanf("%lld %lld %lld %lld %lld", &a, &b, &k, &n, &m) != 5) return 0;

    ll fk = 1, fn = 1, fm = 1;
    for (ll i = 1; i <= k; i++) fk *= i;          // ⚠ k ≥ 21 就溢出了
    for (ll i = 1; i <= n; i++) fn *= i;
    for (ll i = 1; i <= m; i++) fm *= i;
    ll C = fk / (fn * fm);                        // ⚠ 溢出之后这句连整除都不成立

    printf("%lld\n", (((C % P) + P) % P) * qpow(a, n) % P * qpow(b, m) % P);
    return 0;
}

// ★ 第 ③ 版：杨辉三角递推 —— **这一版就已经能 AC 了**
//
// 本章第 3 步那份 brute.cpp 原样搬过来：C(i,j) = C(i−1,j−1) + C(i−1,j)，
// 全程只有加法 ⇒ **不用逆元、不用快速幂、连「p 是不是质数」都不用问**。
//
// k ≤ 1000 ⇒ 三角一共 k(k+1)/2 ≈ 50 万个格子，本机不到 1 毫秒（时限 1 秒）。
// ⇒ ★ 这道题只问**一次**，所以那套 O(k) 预处理 + O(1) 回答的模板在这儿一分钱都不值 ——
//    它是给下一道题（多组询问的 P3807）准备的。
//
// ⚠ 但 a^n、b^m 那两下还是要快速幂（n、m 到 1000，硬乘 1000 次也行，只是没必要）。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll P = 10007;

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

static ll C[1005][1005];

int main() {
    ll a, b, k, n, m;
    if (scanf("%lld %lld %lld %lld %lld", &a, &b, &k, &n, &m) != 5) return 0;

    for (ll i = 0; i <= k; i++) {
        C[i][0] = 1;
        for (ll j = 1; j <= i; j++) C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % P;
    }
    printf("%lld\n", C[k][n] * qpow(a, n) % P * qpow(b, m) % P);
    return 0;
}

// P1313 的度量程序：./p1313Count csv        （本页的数字都出自它）
//
//   line  : ★ 第 ① 版那条溢出线到底在哪 —— 从 k = 0 一路试上去，第一个算错的 k。
//   trig  : ★ 三个错法各自的**触发条件**在四个档里出现几轮（对照「真被抓」看两层差多少）。
//   cost  : ★ 杨辉三角 vs 阶乘+逆元：次数差多少倍、秒表差多少（顶格 k = 1000）。
//
// ⚠ 这里**复刻**了 p1313Gen.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;
const ll P = 10007;

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

/** 复刻 p1313Gen.cpp */
static void gen(unsigned seed, int mode, ll& a, ll& b, ll& k, ll& n, ll& m) {
    mt19937 rng(seed * 1000003u + 20260903u);
    if (mode == 0) { a = (ll)(rng() % 101u); b = (ll)(rng() % 101u); k = (ll)(rng() % 31u); }
    else if (mode == 2) { a = 1; b = 1; k = (ll)(rng() % 1001u); }
    else { a = (ll)(rng() % 1000001u); b = (ll)(rng() % 1000001u); k = (ll)(rng() % 1001u); }
    if (mode == 3) { n = (ll)(rng() % (unsigned)(k + 1)); m = (ll)(rng() % (unsigned)(k + 1)); }
    else { n = (ll)(rng() % (unsigned)(k + 1)); m = k - n; }
}

/** int 版那个快速幂会不会溢出：把它的每一步都拿 long long 走一遍 */
static bool intPowOverflows(ll a, ll e) {
    const ll TOP = 2147483647LL;
    ll res = 1;
    while (e) {
        if (e & 1) { if (res * a > TOP) return true; res = res * a % P; }
        if (a * a > TOP) return true;
        a = a * a % P;
        e >>= 1;
    }
    return false;
}

static ll fac[2005], inv[2005];
static void prep(ll k) {
    fac[0] = 1;
    for (ll i = 1; i <= k; i++) fac[i] = fac[i - 1] * i % P;
    inv[k] = qpow(fac[k], P - 2);
    for (ll i = k; i >= 1; i--) inv[i - 1] = inv[i] * i % P;
}
static ll okAns(ll a, ll b, ll k, ll n, ll m) {
    prep(k);
    ll C = fac[k] * inv[n] % P * inv[k - n] % P;
    return C * qpow(a, n) % P * qpow(b, m) % P;
}
static ll naiveAns(ll a, ll b, ll k, ll n, ll m) {
    ll fk = 1, fn = 1, fm = 1;
    for (ll i = 1; i <= k; i++) fk *= i;
    for (ll i = 1; i <= n; i++) fn *= i;
    for (ll i = 1; i <= m; i++) fm *= i;
    ll C = fk / (fn * fm);
    return (((C % P) + P) % P) * qpow(a, n) % P * qpow(b, m) % P;
}

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

    /* ① 第 ① 版那条线：k 从 0 往上，第一个让它算错的 k（每个 k 把 n 全试一遍） */
    int firstBad = -1;
    for (ll k = 0; k <= 40 && firstBad < 0; k++)
        for (ll n = 0; n <= k; n++)
            if (naiveAns(7, 11, k, n, k - n) != okAns(7, 11, k, n, k - n)) { firstBad = (int)k; break; }

    /* ② 四个档 × 三个触发条件 */
    int trigK[4] = {0, 0, 0, 0}, trigOvf[4] = {0, 0, 0, 0}, trigAB[4] = {0, 0, 0, 0};
    for (int mode = 0; mode < 4; mode++)
        for (int s = 1; s <= 300; s++) {
            ll a, b, k, n, m; gen((unsigned)s, mode, a, b, k, n, m);
            if (k >= 21) trigK[mode]++;                                   // Naive 溢出
            if (intPowOverflows(a, n) || intPowOverflows(b, m)) trigOvf[mode]++;   // NoModA 溢出
            if (a != b) trigAB[mode]++;                                   // Swap 有机会现形
        }

    /* ③ 顶格 k = 1000：杨辉三角 vs 阶乘 + 逆元（次数 ＋ 秒表） */
    const ll K = 1000;
    ll addPascal = K * (K + 1) / 2;          // 三角里真正做加法的格子数
    ll mulFast = K + K + 14;                 // 阶乘 K 次 + 倒推 K 次 + 一次快速幂约 14 步
    static vector<vector<int>> C;
    double msPascal, msFast;
    {
        double t[3];
        for (int r = 0; r < 3; r++) {
            auto t0 = steady_clock::now();
            C.assign(K + 1, vector<int>(K + 1, 0));
            for (ll i = 0; i <= K; i++) { C[i][0] = 1; for (ll j = 1; j <= i; j++) C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % (int)P; }
            t[r] = duration<double, milli>(steady_clock::now() - t0).count();
        }
        msPascal = med3(t[0], t[1], t[2]);
    }
    {
        double t[3];
        for (int r = 0; r < 3; r++) {
            auto t0 = steady_clock::now();
            volatile ll sink = 0;
            for (int rep = 0; rep < 100; rep++) { prep(K); sink += fac[K] + inv[K]; }
            t[r] = duration<double, milli>(steady_clock::now() - t0).count() / 100.0;
        }
        msFast = med3(t[0], t[1], t[2]);
    }

    if (csv) {
        printf("firstBad,%d\n", firstBad);
        for (int i = 0; i < 4; i++) printf("trigK%d,%d\n", i, trigK[i]);
        for (int i = 0; i < 4; i++) printf("trigOvf%d,%d\n", i, trigOvf[i]);
        for (int i = 0; i < 4; i++) printf("trigAB%d,%d\n", i, trigAB[i]);
        printf("addPascal,%lld\n", addPascal);
        printf("mulFast,%lld\n", mulFast);
        printf("opsRatio,%.1f\n", (double)addPascal / (double)mulFast);
        printf("msPascal,%.3f\n", msPascal);
        printf("msFast,%.3f\n", msFast);
        printf("msRatio,%.1f\n", msPascal / msFast);
        printf("pascalUnderLimit,%.0f\n", 1000.0 / msPascal);
        return 0;
    }
    printf("① 第 ① 版（long long 阶乘）那条溢出线\n");
    printf("   k = 0 .. %d 每个 n 都试过 ⇒ 第一个算错的 k = **%d**（因为 21! 装不下 long long）\n", 40, firstBad);
    printf("   ⇒ ★ 题面 30%% 那一档是 k ≤ 10 ⇒ **这一版稳拿 30 分**\n");

    const char* NAME[4] = {"顺手写的（a,b ≤ 100、k ≤ 30）", "照题面（a,b ≤ 10⁶、k ≤ 1000）",
                           "a = b = 1（题面 50% 那一档）", "⚠ 违反题面：n + m ≠ k"};
    printf("\n② 四个档各 300 轮，三个错法的**触发条件**各出现几轮\n");
    for (int i = 0; i < 4; i++)
        printf("   档 %d %-30s k ≥ 21 的 %3d 轮 ／ int 会溢出的 %3d 轮 ／ a ≠ b 的 %3d 轮\n",
               i, NAME[i], trigK[i], trigOvf[i], trigAB[i]);
    printf("   ⇒ ★ 档 2 那两个 0 就是「a = b = 1 ⇒ 弄反了也一样」的证明（题面 50%% 档的形状）\n");

    printf("\n③ 顶格 k = 1000：两种正确写法（独占实测，3 次取中位数；时限 1000 ms）\n");
    printf("   杨辉三角        %8lld 次加法    %.3f ms   ⇒ 余量 %.0f 倍\n", addPascal, msPascal, 1000.0 / msPascal);
    printf("   阶乘 + 逆元     %8lld 次乘法    %.3f ms   ⇒ 次数少 %.1f 倍、秒表快 %.1f 倍\n",
           mulFast, msFast, (double)addPascal / (double)mulFast, msPascal / msFast);
    printf("   ⇒ ★ 两版都远在时限之内 ⇒ **这道题该交的是杨辉三角**，逆元那套是给下一道题练手的\n");
    return 0;
}