37 · 最大公约数、同余与快速幂

上一节把数拆成质数，这一节处理“两个数的公共结构”和“大指数”。三件事：两个数的最大公约数与最小公倍数、同余方程与逆元、a^b mod m。第三件是第二件的工具，第二件是第一件的推广——三者用同一套代码互相支撑。

先把这几个工具拆开看清。它们都很短，但每一个都有一处“不这么写就错”的细节，这些细节正是这一节真正的考点。

#include <bits/stdc++.h>
using namespace std;

int main() {
    cout << gcd(12, 18) << "\n";        // 6
    cout << gcd(0, 7)   << "\n";        // 7：gcd(0, x) = |x|
    cout << gcd(-12, 18)<< "\n";        // 6：取绝对值

    long long a = 1000000000000LL, b = 8;
    cout << gcd(a, b) << "\n";          // 4：两个都是 long long

    // 手写辗转相除（不依赖 C++17）
    long long x = a, y = b;
    while (y) { long long t = x % y; x = y; y = t; }
    cout << x << "\n";                  // 4
    return 0;
}

long long lcm(long long a, long long b) {
    return a / gcd(a, b) * b;      // 先除后乘
}

// 反例：a = b = 4e9
//   a * b / gcd = 1.6e19 / 4e9   → 1.6e19 已超出 long long，溢出成负数
//   a / gcd * b = 1e9 * 4e9      → 4e18，正确

// ① 防溢出的乘法取模：a、b < m ≤ 10^18
long long mulMod(long long a, long long b, long long m) {
    return (long long)((__int128)a * b % m);   // 中间 10^36，long long 装不下
}

// ② 打印 __int128：cout 不认识它，只能自己按位输出
void print128(__int128 x) {
    if (x < 0) { putchar('-'); x = -x; }
    if (x > 9) print128(x / 10);
    putchar('0' + (int)(x % 10));
}

long long qpow(long long a, long long b, long long mod) {
    long long r = 1 % mod;              // 处理 mod = 1 的退化情况
    a %= mod;
    while (b) {
        if (b & 1) r = r * a % mod;     // 这一位是 1：把 a^(2^k) 乘进来
        a = a * a % mod;                // a 平方，准备下一位
        b >>= 1;
    }
    return r;
}

// qpow(2, 10, 9)   = 1024 % 9 = 7
// qpow(2, 0, 1000) = 1
// qpow(7, 5, 1)    = 0

r = 1, a = b % k
当 p > 0：
    若 p 的最低位是 1：r = r * a % k     ← 这一位对应的幂要乘进来
    a = a * a % k                        ← 幂次翻倍：a, a², a⁴, ...
    p >>= 1                              ← 右移，看下一位
输出 r

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

// 计算 a^b mod m
ll qpow(ll a, ll b, ll m) {
    ll r = 1 % m;                          // m = 1 时答案恒为 0，1 % 1 = 0
    a %= m;                                // 底数先进模，后面的乘法才安全
    while (b) {
        if (b & 1) r = (__int128)r * a % m; // 这一位是 1，把当前的 a^(2^k) 乘进答案
        a = (__int128)a * a % m;            // 平方，得到下一份的幂
        b >>= 1;                            // 处理下一个二进制位
    }
    return r;
}

int main() {
    ll b, p, k;
    cin >> b >> p >> k;                     // p 可达 2^31，必须 long long
    cout << b << "^" << p << " mod " << k << "=" << qpow(b, p, k) << "\n";
    return 0;
}

快速幂解决“指数很大怎么算”。数论的另一半问题是“除法怎么取模”。加减乘都可以边算边取模，只有除法不行：(a / b) mod m 既不等于 (a mod m) / (b mod m)，很多时候压根不是整数。要让“除以 b”变成“乘以 b 的逆元”，就得解同余方程 b·x ≡ 1 (mod m)。

求逆元要先会解 a·x + b·y = gcd(a, b)。裴蜀定理保证这个方程一定有整数解，而把辗转相除的每一步“回代”就能把解构造出来——这就是扩展欧几里得。

// 求 a*x + b*y = gcd(a,b) 的一组整数解，返回 gcd(a,b)
// x、y 用引用传递，返回值通过参数带出（第 11 节学过）
long long exgcd(long long a, long long b, long long &x, long long &y) {
    if (b == 0) {                    // 出口：a*1 + 0*0 = a
        x = 1; y = 0;
        return a;
    }
    long long x1, y1;
    long long g = exgcd(b, a % b, x1, y1);   // 先让规模变小，再回代
    x = y1;                          // a 的系数
    y = x1 - (a / b) * y1;           // b 的系数
    return g;
}

拿到 exgcd 之后，逆元就是一行的事：a·x + m·y = 1 里的 x 就是 a 模 m 的逆元。另一种求法是把上一节的快速幂接上来——当 m 是质数时，费马小定理直接给答案。

// ① 扩展欧几里得：任意模数，只要 gcd(a, m) == 1 就有逆元
long long invExgcd(long long a, long long m) {
    long long x, y;
    long long g = exgcd(a, m, x, y);     // a*x + m*y = gcd(a, m)
    if (g != 1) return -1;               // 不互质，逆元不存在
    return (x % m + m) % m;              // 归一化到 [0, m)
}

// ② 费马小定理：只在 m 是质数时可用（1e9+7、998244353 这类）
long long invFermat(long long a, long long m) {
    return qpow(a, m - 2, m);            // a^(m-2) ≡ a^(-1) (mod m)
}

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

// 求 a*x + b*y = gcd(a,b) 的一组整数解，返回 gcd(a,b)
ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) {                       // 出口：a*1 + 0*0 = a
        x = 1; y = 0;
        return a;
    }
    ll x1, y1;
    ll g = exgcd(b, a % b, x1, y1);     // 规模先变小
    x = y1;                             // 回代：a 的系数
    y = x1 - (a / b) * y1;              // 回代：b 的系数
    return g;
}

int main() {
    ll a, b, x, y;
    cin >> a >> b;

    exgcd(a, b, x, y);                  // 得到 a*x + b*y = gcd(a,b) = 1
    x = (x % b + b) % b;                // 归一化到 [0, b)：先 % 再 +b 再 %
    cout << x << "\n";
    return 0;
}