洛谷 P1082 [NOIP 2012 提高组] 同余方程

// 第 ① 版：从 1 试到 b−1，看谁满足 a·x ≡ 1 —— 参照物，而且它值 60 分
//
// ★ 题面三档：
//   40%：b ≤ 10³      ⇒ 一千次，秒过；
//   60%：b ≤ 5×10⁷    ⇒ 五千万次乘法取模，本机几百毫秒 ⇒ **它拿得到 60 分**；
//   100%：b ≤ 2×10⁹   ⇒ 二十亿次 ⇒ 超时。
// ⇒ 又一次「题面的分档是出题人递过来的工具」。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

int main() {
    ll a, b;
    if (scanf("%lld %lld", &a, &b) != 2) return 0;
    a %= b;
    for (ll x = 1; x < b; x++)
        if (a * x % b == 1) { printf("%lld\n", x); return 0; }
    printf("-1\n");            // 题面保证有解，走不到这里
    return 0;
}

// P1082 同余方程 —— 正解：扩展欧几里得求逆元
//
// ★★ 题面这一句是全部的分水岭：**b 不保证是质数。**
//   · b 是质数 ⇒ 费马小定理，逆元 = a^(b−2) mod b，一次快速幂就完了；
//   · b 不保证 ⇒ **费马用不了**，只能 exgcd（[第 40 章第 13 步](/ch/40-gcd-euclid/)那份）。
//
// ax ≡ 1 (mod b) ⇔ 存在 y 使 a·x + b·y = 1。
// 而 exgcd 解的正是 a·x + b·y = gcd(a, b) —— 题面保证有解 ⇒ gcd(a, b) = 1 ⇒ 直接用。
//
// ⚠ 最后一步不能省：exgcd 给的 x 可能是负数（也可能大于 b），
//   而题目要的是**最小正整数解** ⇒ ((x % b) + b) % b。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

/** 返回 gcd(a,b)，并让 a·x + b·y = gcd */
static ll exgcd(ll a, ll b, ll& x, ll& y) {
    if (!b) { x = 1; y = 0; return a; }
    ll x1, y1;
    ll g = exgcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - a / b * y1;
    return g;
}

int main() {
    ll a, b, x, y;
    if (scanf("%lld %lld", &a, &b) != 2) return 0;
    exgcd(a, b, x, y);
    printf("%lld\n", (x % b + b) % b);
    return 0;
}

// P1082 的度量程序：./p1082Count csv        （本页的数字都出自它）
//
//   int   : ★★★ 全用 int 的那一版**到底哪一步溢出** —— 把它的每一步都拿 long long 复算一遍。
//           ⚠ 结论和直觉不一样：**exgcd 的中间系数是安全的，挂的是最后那句归正。**
//   coef  : ★ 上一条的证据：exgcd 一路回代出来的系数，绝对值最大是多少（拿 max(a,b) 比一比）。
//   prime : ⚠ 随机挑一个 b，它是质数的概率 —— 这就是「费马版被顺手写的生成器盖住」的量。
//   verify: ★ 顶格那一档没法用暴力当参照物 ⇒ 改用**验证器**：a·x ≡ 1 (mod b) 且 1 ≤ x < b。
//           （满足这两条的 x 唯一，所以验证器和对拍一样强，而且便宜得多。）
//   ms    : 题面三档 × 两种写法的秒表。
//
// ⚠ 这里**复刻**了 p1082Gen.cpp 的档位逻辑（同一个 mt19937、同一个种子公式）。
// ⚠ 秒表用 steady_clock 在进程内量；模数取运行期的值（[P1965 那一跤](/sol/p1965/)：
//   编译期常量会让 -O2 把取模换成乘法）。
#include <bits/stdc++.h>
#include <chrono>
using namespace std;
using namespace std::chrono;
typedef long long ll;

static const ll INT_TOP = 2147483647LL;

static ll gcdll(ll a, ll b) { while (b) { ll t = a % b; a = b; b = t; } return a; }
static bool isPrime(ll x) { if (x < 2) return false; for (ll i = 2; i * i <= x; i++) if (x % i == 0) return false; return true; }

/** 标准 exgcd（long long），顺手记下这一路系数的最大绝对值 */
static ll maxCoef;
static ll exgcd(ll a, ll b, ll& x, ll& y) {
    if (!b) { x = 1; y = 0; return a; }
    ll x1, y1;
    ll g = exgcd(b, a % b, x1, y1);
    x = y1; y = x1 - a / b * y1;
    maxCoef = max(maxCoef, max(llabs(x), llabs(y)));
    return g;
}

/** 把 int 版的每一步都用 long long 复算：分别看 exgcd 内部和最后那句会不会越界 */
static bool exgcdIntBad(ll a, ll b, ll& x, ll& y) {
    if (!b) { x = 1; y = 0; return false; }
    ll x1, y1;
    bool bad = exgcdIntBad(b, a % b, x1, y1);
    ll q = a / b;
    if (llabs(q * y1) > INT_TOP) bad = true;
    ll ny = x1 - q * y1;
    if (llabs(ny) > INT_TOP) bad = true;
    x = y1; y = ny;
    return bad;
}

static void gen(unsigned seed, int mode, ll& a, ll& b) {
    mt19937 rng(seed * 1000003u + 20260902u);
    while (true) {
        if (mode == 0) { a = 2 + (ll)(rng() % 999u); b = 2 + (ll)(rng() % 999u); }
        else if (mode == 1) { b = 2 + (ll)(rng() % 999u); if (!isPrime(b)) continue; a = 2 + (ll)(rng() % 999u); }
        else if (mode == 2) { b = 4 + (ll)(rng() % 997u); if (isPrime(b)) continue; a = 2 + (ll)(rng() % 999u); }
        else { a = 1000000000LL + (ll)(rng() % 1000000000u); b = 1000000000LL + (ll)(rng() % 1000000000u); }
        if (b >= 2 && gcdll(a, b) == 1) break;
    }
}
static double med3(double a, double b, double c) { return max(min(a, b), min(max(a, b), c)); }

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

    /* ① 四个档：int 在 exgcd 内部越界的轮数 / 在最后那句越界的轮数 / x 是负数的轮数 / 验证器 */
    int badIn[4] = {0, 0, 0, 0}, badTail[4] = {0, 0, 0, 0}, negX[4] = {0, 0, 0, 0}, verifyBad[4] = {0, 0, 0, 0};
    ll coefOverMax[4] = {0, 0, 0, 0};
    for (int mode = 0; mode < 4; mode++)
        for (int s = 1; s <= 300; s++) {
            ll a, b, x, y; gen((unsigned)s, mode, a, b);
            maxCoef = 0;
            exgcd(a, b, x, y);
            if (maxCoef > coefOverMax[mode]) coefOverMax[mode] = maxCoef;
            if (x < 1) negX[mode]++;
            ll xi, yi;
            if (exgcdIntBad(a, b, xi, yi)) badIn[mode]++;
            ll r = xi % b;                       // int 版最后那句：(x % b + b) % b
            if (llabs(r) > INT_TOP || llabs(r + b) > INT_TOP) badTail[mode]++;
            ll ans = (x % b + b) % b;            // 正解
            if (!(ans >= 1 && ans < b && (__int128)a % b * ans % b == 1)) verifyBad[mode]++;
        }

    /* ② exgcd 的系数到底有多大：拿 max(a,b) 当尺子，看比值 */
    double coefRatio = 0;
    {
        mt19937 rng(20260902u);
        for (int i = 0; i < 200000; i++) {
            ll a = 1 + (ll)(rng() % 2000000000u), b = 1 + (ll)(rng() % 2000000000u);
            if (b < 2 || gcdll(a, b) != 1) continue;
            ll x, y; maxCoef = 0; exgcd(a, b, x, y);
            coefRatio = max(coefRatio, (double)maxCoef / (double)max(a, b));
        }
    }

    /* ③ 随机 b ≤ 1000 里质数占多少 */
    int primeCnt = 0;
    for (int b = 2; b <= 1000; b++) if (isPrime(b)) primeCnt++;

    /* ④ 秒表：暴力在三档顶格 ／ 正解 */
    // ⚠⚠ 暴力是「找到就 break」的 ⇒ **必须让它扫满整段才算最坏**。
    //   第一版顺手写 a = 999999937 % b，那个 a 的逆元恰好很小，量出来 53 毫秒 ——
    //   量的是「运气好的那一次」。这里改成：先挑一个很大的目标 x₀ = b − 2，
    //   再反解出以它为逆元的 a（a ≡ x₀⁻¹），于是循环一定要走到 b − 2。
    static volatile ll VB1 = 1000, VB2 = 50000000LL;
    double msB40, msB60, projB100, usFast;
    auto worstA = [](ll b) {                     // 返回一个 a，使 a 的逆元 = b − 2
        ll x, y; exgcd(b - 2, b, x, y);
        return (x % b + b) % b;
    };
    {
        double t[3];
        for (int r = 0; r < 3; r++) {
            ll b = VB1, a = worstA(b);
            auto t0 = steady_clock::now();
            volatile ll got = 0;
            for (ll x = 1; x < b; x++) if (a * x % b == 1) { got = x; break; }
            (void)got;
            t[r] = duration<double, milli>(steady_clock::now() - t0).count();
        }
        msB40 = med3(t[0], t[1], t[2]);
    }
    {
        ll b = VB2, a = worstA(b);               // 让它扫满整段
        auto t0 = steady_clock::now();
        volatile ll got = 0;
        for (ll x = 1; x < b; x++) if (a * x % b == 1) { got = x; break; }
        (void)got;
        double ms = duration<double, milli>(steady_clock::now() - t0).count();
        msB60 = ms;
        projB100 = ms * (2000000000.0 / 50000000.0);
    }
    {
        auto t0 = steady_clock::now();
        // ⚠ a 和 b 都得随循环变，否则 -O2 会把整次调用提到循环外（[P1965 那一跤](/sol/p1965/)）
        static volatile ll VA0 = 1999999973LL, VB0 = 1999999943LL;
        volatile ll acc = 0;
        for (int r = 0; r < 20000; r++) {
            ll a = VA0 - r, b = VB0 - (r & 1), x, y;
            exgcd(a, b, x, y);
            acc = acc + ((x % b + b) % b);
        }
        usFast = duration<double, micro>(steady_clock::now() - t0).count() / 20000.0;
    }

    if (csv) {
        for (int i = 0; i < 4; i++) printf("badIn%d,%d\n", i, badIn[i]);
        for (int i = 0; i < 4; i++) printf("badTail%d,%d\n", i, badTail[i]);
        for (int i = 0; i < 4; i++) printf("negX%d,%d\n", i, negX[i]);
        for (int i = 0; i < 4; i++) printf("verifyBad%d,%d\n", i, verifyBad[i]);
        printf("coefRatio,%.3f\n", coefRatio);
        printf("primeCnt,%d\n", primeCnt);
        printf("primePct,%.1f\n", primeCnt * 100.0 / 999.0);
        printf("msB40,%.2f\n", msB40);
        printf("msB60,%.0f\n", msB60);
        printf("projB100,%.0f\n", projB100);
        printf("usFast,%.3f\n", usFast);
        return 0;
    }
    const char* NAME[4] = {"顺手写的（≤1000）", "b 取质数", "b 取合数", "顶格（10⁹~2×10⁹）"};
    printf("① 四个档各 300 轮\n");
    for (int i = 0; i < 4; i++)
        printf("   档 %d %-20s int 在 exgcd 内部越界 %3d 轮 ／ **在最后那句归正越界 %3d 轮** ／ x < 1 的 %3d 轮 ／ 验证器不通过 %d 轮\n",
               i, NAME[i], badIn[i], badTail[i], negX[i], verifyBad[i]);
    printf("   ⇒ ★★★ int 挂的地方**不在 exgcd 里**，在 `(x %% b + b) %% b` 那一句：\n");
    printf("      x %% b 最大可到 b−1 ≈ 2×10⁹，再加一个 b 就是 4×10⁹ —— 而 int 只到 2.1×10⁹\n");
    printf("\n② exgcd 的系数有多大：随机二十万组里，|系数| / max(a,b) 的最大值 = %.3f\n", coefRatio);
    printf("   ⇒ 系数被 max(a, b) 界住 ⇒ 题面 a, b ≤ 2×10⁹ 时它们**本身**是塞得进 int 的（余量 7.4%%）\n");
    printf("\n③ 2 ≤ b ≤ 1000 里质数有 %d 个（%.1f%%）—— 顺手挑一个 b，它是质数的概率\n",
           primeCnt, primeCnt * 100.0 / 999.0);
    printf("\n④ 秒表（独占实测；时限 1000 ms）\n");
    printf("   暴力 b = 10³（40%% 档顶格，扫满）  %.2f ms ⇒ 稳拿 40 分\n", msB40);
    printf("   暴力 b = 5×10⁷（60%% 档顶格，扫满）%.0f ms ⇒ 稳拿 60 分\n", msB60);
    printf("   暴力 b = 2×10⁹（100%% 档，外推）   %.0f ms ⇒ ✗ 超时\n", projB100);
    printf("   正解 exgcd（顶格，两万次取平均）   %.2f 微秒 ⇒ ★ AC\n", usFast);
    return 0;
}