最大公约数（GCD）：从原理到实战

在算法竞赛中，数论是不可或缺的一部分，而最大公约数（Greatest Common Divisor, GCD）则是数论问题的基石。无论是直接考察，还是作为解决更复杂问题的工具，深刻理解 GCD 的原理、高效计算方法及其应用都至关重要。本文将带你深入探索 GCD 的世界。

1. 核心概念

1.1 定义

对于两个非零整数 a 和 b，它们的最大公约数 gcd(a, b) 是能够同时整除 a 和 b 的最大正整数。

形式化定义： d = gcd(a, b) 当且仅当

1. d | a 且 d | b （d 是 a 和 b 的公约数）
2. 对于任意满足 c | a 且 c | b 的整数 c，都有 c | d （d 是所有公约数中最大的）

特殊约定：gcd(a, 0) = |a|。

1.2 基本性质

• gcd(a, b) = gcd(b, a) （交换律）
• gcd(a, b) = gcd(|a|, |b|) （我们通常处理正整数）
• gcd(a, a) = |a|
• gcd(a, 0) = |a|
• gcd(a, ka) = |a| （对于任意整数 k）
• gcd(a, b) = gcd(a - b, b) （更相减损术的核心，也是欧几里得算法的基石）
• gcd(a, b, c) = gcd(gcd(a, b), c) （结合律，可推广到多个数）
• gcd(ka, kb) = k * gcd(a, b) （k > 0）
• 若 gcd(a, b) = d，则 gcd(a/d, b/d) = 1。称 a/d 和 b/d 互质。

理解这些性质是灵活运用 GCD 的前提。

2. 核心算法：欧几里得算法 (Euclidean Algorithm)

计算 GCD 最常用且高效的方法是欧几里得算法，也称为辗转相除法。

2.1 算法原理

该算法基于一个关键定理： 对于任意两个非负整数 a 和 b（不全为零），有：

$$gcd(a, b) = gcd(b, a \pmod b)$$

其中 a mod b 表示 a 除以 b 的余数。

证明思路: 设 d = gcd(a, b)。则 d | a 且 d | b。 令 a = kb + r，其中 r = a mod b，0 <= r < b。 因为 d | a 且 d | b，所以 d | (a - kb)，即 d | r。 这说明 d 也是 b 和 r 的公约数。

设 d' = gcd(b, r)。则 d' | b 且 d' | r。 因为 a = kb + r，所以 d' | (kb + r)，即 d' | a。 这说明 d' 也是 a 和 b 的公约数。

因为 d 是 a 和 b 的最大公约数，且 d' 是 a 和 b 的公约数，所以 d' <= d。 因为 d' 是 b 和 r 的最大公约数，且 d 是 b 和 r 的公约数，所以 d <= d'。 综上，d = d'，即 gcd(a, b) = gcd(b, a mod b)。

2.2 递归实现

基于上述原理，可以轻松写出递归代码。

#include <bits/stdc++.h>

using namespace std;
using ll = long long;

// 递归实现欧几里得算法
// 计算 a 和 b 的最大公约数
// 时间复杂度: O(log(min(a, b)))
ll gcd(ll a, ll b) {
    // 基本情况：当 b 为 0 时，gcd(a, 0) = a
    // 递归终止条件
    return b == 0 ? a : gcd(b, a % b);
}

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

    ll a, b;
    cin >> a >> b;
    cout << gcd(a, b) << endl;

    return 0;
}

代码解释:

• 函数 gcd(a, b) 计算 a 和 b 的最大公约数。
• 递归的终止条件是 b == 0，此时 a 就是最大公约数。
• 递归步骤是 gcd(b, a % b)，不断用除数和余数替换原来的两个数。
• 由于每次取模操作至少会使其中一个数减半（粗略估计），因此递归深度是对数级别的。

复杂度分析:

• 时间复杂度: O(log(min(a, b)))。这是一个非常高效的算法。最坏情况发生在输入为连续的斐波那契数时，但仍然是对数复杂度。
• 空间复杂度: O(log(min(a, b))) （递归调用栈的深度）。

2.3 迭代实现

在竞赛中，有时为了避免递归深度过大或追求极致的常数效率，会使用迭代实现。

#include <bits/stdc++.h>

using namespace std;
using ll = long long;

// 迭代实现欧几里得算法
// 计算 a 和 b 的最大公约数
// 时间复杂度: O(log(min(a, b)))
ll gcd(ll a, ll b) {
    while (b != 0) {
        ll r = a % b;
        a = b;
        b = r;
    }
    return a; // 循环结束时，a 保存着 gcd
}

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

    ll a, b;
    cin >> a >> b;
    cout << gcd(a, b) << endl;

    return 0;
}

代码解释:

• 使用 while 循环代替递归。
• 循环条件是 b != 0。
• 在循环内部，计算余数 r，然后更新 a 和 b，相当于 (a, b) -> (b, a % b) 的转换。
• 当 b 变为 0 时，循环终止，此时的 a 即为所求 GCD。

复杂度分析:

• 时间复杂度: O(log(min(a, b)))。
• 空间复杂度: O(1)。迭代版本空间效率更高。

在 C++ STL numeric 头文件中，也提供了 std::gcd 函数，可以直接使用。

3. 扩展欧几里得算法 (Extended Euclidean Algorithm - exGCD)

扩展欧几里得算法不仅计算 gcd(a, b)，还能找到一对整数 x 和 y，使得它们满足贝祖定理（Bézout's identity）：

$$ax + by = gcd(a, b)$$

3.1 算法原理

exGCD 同样基于 gcd(a, b) = gcd(b, a mod b)。 假设我们已经求得了一组 x' 和 y' 满足：

$$bx' + (a \pmod b)y' = gcd(b, a \pmod b)$$

因为 gcd(a, b) = gcd(b, a mod b)，所以：

$$bx' + (a \pmod b)y' = gcd(a, b)$$

我们知道 a mod b = a - floor(a / b) * b。代入上式：

$$bx' + (a - \lfloor a/b \rfloor * b)y' = gcd(a, b)$$

整理得到：

$$ay' + b(x' - \lfloor a/b \rfloor y') = gcd(a, b)$$

对比 ax + by = gcd(a, b)，我们可以得到：

$$x = y'$$ $$y = x' - \lfloor a/b \rfloor y'$$

这给出了从 (x', y') 推导到 (x, y) 的递推关系。

递归的边界条件是当 b = 0 时，gcd(a, 0) = a。此时方程变为 ax + 0y = a。显然，一组特解是 x = 1, y = 0。

3.2 实现

#include <bits/stdc++.h>

using namespace std;
using ll = long long;

// 扩展欧几里得算法
// 计算 ax + by = gcd(a, b) 的一组解 (x, y)
// 返回值为 gcd(a, b)
// x 和 y 是通过引用传递的参数，用于接收结果
// 时间复杂度: O(log(min(a, b)))
ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a; // 边界情况 gcd(a, 0) = a, 此时 ax + 0y = a 的解为 x=1, y=0
    }
    ll x1, y1;
    // 递归调用，求解 b*x1 + (a%b)*y1 = gcd(b, a%b)
    ll d = exgcd(b, a % b, x1, y1);
    // 根据推导关系计算当前的 x, y
    x = y1;
    y = x1 - (a / b) * y1;
    return d; // 返回 gcd(a, b)
}

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

    ll a, b, x, y;
    cin >> a >> b;

    ll d = exgcd(a, b, x, y);

    cout << "gcd(" << a << ", " << b << ") = " << d << endl;
    cout << "x = " << x << ", y = " << y << endl;
    cout << a << "*" << x << " + " << b << "*" << y << " = " << a * x + b * y << endl; // 验证

    return 0;
}

代码解释:

• 函数 exgcd(a, b, x, y) 通过引用参数 x 和 y 返回贝祖等式的一组解，并返回 gcd(a, b)。
• 递归基是 b == 0，此时 x=1, y=0，返回 a。
• 递归调用 exgcd(b, a % b, x1, y1) 得到下一层的解 x1, y1 和 gcd 值 d。
• 根据推导公式 x = y1 和 y = x1 - (a / b) * y1 计算当前层的解 x, y。
• 注意 a / b 在 C++ 中对于整数是整除，即 floor(a / b)。

复杂度分析:

• 时间复杂度: O(log(min(a, b)))，与普通欧几里得算法相同。
• 空间复杂度: O(log(min(a, b))) (递归栈)。同样可以写成迭代形式以优化空间，但代码会稍显复杂。

3.3 应用：求解线性丢番图方程

形如 ax + by = c 的方程称为线性丢番图方程，其中 a, b, c 是整数，要求解整数 x, y。

定理: 线性丢番图方程 ax + by = c 有整数解 (x, y) 当且仅当 gcd(a, b) 整除 c。

求解步骤:

1. 计算 d = gcd(a, b)。
2. 检查 c % d == 0。如果不是，方程无整数解。
3. 如果 c % d == 0，则先用 exgcd(a, b, x0, y0) 求出 ax0 + by0 = d 的一组特解 (x0, y0)。
4. 将等式两边同乘以 c / d： a(x0 * c / d) + b(y0 * c / d) = d * c / d = c
5. 于是，原方程 ax + by = c 的一组特解是： x = x0 * (c / d) y = y0 * (c / d)
6. 通解: 如果 (x, y) 是一组特解，则所有整数解可以表示为： x_k = x + k * (b / d) y_k = y - k * (a / d) 其中 k 是任意整数。

代码示例框架:

// 假设 exgcd 函数已定义
// 求解 ax + by = c
void solve_diophantine(ll a, ll b, ll c) {
    ll x0, y0;
    ll d = exgcd(a, b, x0, y0);

    if (c % d != 0) {
        cout << "No integer solution\n";
    } else {
        ll k = c / d;
        ll x = x0 * k;
        ll y = y0 * k;
        // x, y 是 ax + by = c 的一组特解

        // 如果需要求最小正整数解 x 或其他特定解，
        // 可以利用通解公式调整
        // 例如求最小非负整数解 x:
        ll b_d = b / d;
        x = (x % b_d + b_d) % b_d; // C++负数取模特性处理
        // 相应地调整 y (需要根据具体题目要求)
        // y = (c - a * x) / b; // 确保 b 不为 0

        cout << "One solution: x = " << x << ", y = " << y << "\n";
        // 通解形式: x_k = x + k*(b/d), y_k = y - k*(a/d)
    }
}

3.4 应用：求解模线性方程和乘法逆元

求解形如 ax ≡ c (mod m) 的模线性方程。 这等价于求解 ax + my = c 的整数解 (x, y)，其中 y 是某个整数。 这又回到了线性丢番图方程 ax + my = c。

步骤:

1. 令 d = gcd(a, m)。
2. 如果 c 不能被 d 整除，即 c % d != 0，则方程 ax ≡ c (mod m) 无解。
3. 如果 c % d == 0，则方程有 d 个模 m 意义下不同的解。
   • 先用 exgcd(a, m, x0, y0) 求得 ax0 + my0 = d 的解 (x0, y0)。
   • 原方程 ax ≡ c (mod m) 的一个解可以通过 x = x0 * (c / d) mod m 得到。但要注意，这只是其中一个解，且需要考虑负数取模。
   • 更稳妥的方式是：将方程两边同除以 d，得到 (a/d)x ≡ (c/d) (mod m/d)。
   • 令 a' = a/d, c' = c/d, m' = m/d。此时 gcd(a', m') = 1。
   • 问题转化为求解 a'x ≡ c' (mod m')。
   • 因为 gcd(a', m') = 1，所以 a' 模 m' 存在乘法逆元。
   • 用 exgcd(a', m', inv, y) 求解 a' * inv + m' * y = 1，得到 inv，即 a' 模 m' 的逆元。
   • 两边同乘以 inv：inv * a' * x ≡ inv * c' (mod m')
   • 得到 x ≡ inv * c' (mod m')。
   • 设 x_ans = (inv * c') % m'。为了得到正数解，可以写成 x_ans = (inv * c' % m' + m') % m'。
   • 这个 x_ans 是模 m' 意义下的唯一解。
   • 原方程 ax ≡ c (mod m) 的所有 d 个解为： x_k = x_ans + k * m'，其中 k = 0, 1, ..., d-1。

乘法逆元: 当我们需要求解 ax ≡ 1 (mod m) 时，就相当于求 a 模 m 的乘法逆元 x。 这要求 ax + my = 1 有解。根据线性丢番图方程的理论，这当且仅当 gcd(a, m) = 1。 此时，直接调用 exgcd(a, m, x, y)，得到的解 x 就是 a 模 m 的一个逆元。 为了得到 [0, m-1] 范围内的最小正整数逆元，通常取 (x % m + m) % m。

代码示例：求逆元

#include <bits/stdc++.h>

using namespace std;
using ll = long long;

// 扩展欧几里得算法 (同上)
ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) {
        x = 1; y = 0; return a;
    }
    ll x1, y1, d = exgcd(b, a % b, x1, y1);
    x = y1; y = x1 - (a / b) * y1;
    return d;
}

// 计算 a 模 m 的乘法逆元
// 前提：gcd(a, m) = 1
// 返回值: 若逆元存在，返回最小正整数逆元；否则返回 -1 (或抛出异常)
ll modInverse(ll a, ll m) {
    ll x, y;
    ll d = exgcd(a, m, x, y);
    if (d != 1) {
        return -1; // 逆元不存在
    }
    // 确保返回的是 [0, m-1] 范围内的最小正整数
    return (x % m + m) % m;
}

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

    ll a, m;
    cout << "Enter a and m to find modular inverse of a mod m:\n";
    cin >> a >> m;

    ll inv = modInverse(a, m);

    if (inv == -1) {
        cout << "Inverse does not exist (gcd(a, m) != 1)\n";
    } else {
        cout << "Modular inverse of " << a << " mod " << m << " is " << inv << endl;
        cout << "Verification: (" << a << " * " << inv << ") % " << m << " = " << (a * inv) % m << endl;
    }

    return 0;
}

4. 最小公倍数 (Least Common Multiple - LCM)

最小公倍数 lcm(a, b) 是能够被 a 和 b 同时整除的最小正整数。 GCD 和 LCM 之间有一个重要的关系：

$$gcd(a, b) \times lcm(a, b) = |a \times b|$$

因此，可以通过 GCD 来计算 LCM：

$$lcm(a, b) = \frac{|a \times b|}{gcd(a, b)}$$

注意潜在的溢出问题： 计算 a * b 时可能会超出整数类型的表示范围。为了避免这种情况，可以使用更安全的计算方式：

$$lcm(a, b) = \frac{|a|}{gcd(a, b)} \times |b|$$

或者

$$lcm(a, b) = |a| \times \frac{|b|}{gcd(a, b)}$$

这样可以先做除法，降低中间结果的大小。

代码实现:

#include <bits/stdc++.h>

using namespace std;
using ll = long long;

// 欧几里得算法 (同上)
ll gcd(ll a, ll b) {
    return b == 0 ? a : gcd(b, a % b);
}

// 计算 a 和 b 的最小公倍数
// 使用安全的计算方式避免溢出
// 时间复杂度: O(log(min(a, b))) (主要来自 gcd 计算)
ll lcm(ll a, ll b) {
    if (a == 0 || b == 0) return 0; // 处理边界情况
    // a / gcd(a, b) * b
    // 确保使用 64 位整数 (long long)
    return (a / gcd(a, b)) * b;
}

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

    ll a, b;
    cin >> a >> b;
    // 处理负数输入，通常 lcm 是针对正整数的
    ll abs_a = abs(a);
    ll abs_b = abs(b);

    cout << "lcm(" << abs_a << ", " << abs_b << ") = " << lcm(abs_a, abs_b) << endl;

    return 0;
}

5. GCD 的其他性质与应用

5.1 多个数的 GCD 和 LCM

• GCD: 利用结合律 gcd(a, b, c) = gcd(gcd(a, b), c)。可以逐个计算： result = gcd(a1, a2) result = gcd(result, a3) ... result = gcd(result, an)

  // 计算数组 nums 中所有元素的 GCD
  ll multi_gcd(const vector<ll>& nums) {
      if (nums.empty()) return 0; // 或者根据题目定义返回
      ll res = nums[0];
      for (size_t i = 1; i < nums.size(); ++i) {
          res = gcd(res, nums[i]);
          if (res == 1) break; // 优化：如果 gcd 已经是 1，后续无需计算
      }
      return res;
  }
  

• LCM: lcm(a, b, c) 不能简单地用 lcm(lcm(a, b), c) 来套用 a*b/gcd 公式，因为会涉及重复除以公因子。 正确的关系通常与素数分解相关。 但在实践中，可以利用 lcm(a,b) = a/gcd(a,b)*b 迭代计算： result = lcm(a1, a2) result = lcm(result, a3) ... result = lcm(result, an) 注意：迭代计算 LCM 时，中间结果可能会非常大，极易溢出，需要特别注意数据范围或使用高精度计算。

  // 计算数组 nums 中所有元素的 LCM
  // 极易溢出！仅适用于结果在 long long 范围内的情况
  ll multi_lcm(const vector<ll>& nums) {
      if (nums.empty()) return 0;
      ll res = nums[0];
      for (size_t i = 1; i < nums.size(); ++i) {
          // 需要检查中间结果是否会溢出
          // 这里简化处理，假设不会溢出或溢出行为符合预期（例如自然溢出）
          res = lcm(res, nums[i]);
      }
      return res;
  }
  

5.2 类欧几里得算法 (Generalized Euclidean Algorithm)

类欧几里得算法用于计算形如下式的和：

$$\sum_{i=0}^{n} \lfloor \frac{ai + b}{c} \rfloor$$

以及类似的带有 floor 的求和式子。其推导过程与欧几里得算法类似，通过递归将问题规模缩小。这是一个相对高级的内容，通常出现在较难的数论题目中。此处不展开细节，但需要知道有这类算法与 GCD 的思想相关。

5.3 GCD 与数据结构结合：区间 GCD

有时需要快速查询一个区间内所有数的最大公约数。GCD 满足结合律且具有可合并性，因此可以使用线段树或稀疏表（ST表）来维护区间 GCD。

线段树维护区间 GCD:

• 节点信息: 每个节点存储对应区间的 GCD 值。
• PushUp 操作: tree[p].gcd = gcd(tree[p<<1].gcd, tree[p<<1|1].gcd)。父节点的 GCD 等于左右子节点 GCD 的 GCD。
• Build 操作: 叶子节点存储单个元素值，然后自底向上 pushup。
• Query 操作: 标准线段树区间查询，合并查询到的子区间 GCD 结果。

代码框架 (线段树):

#include <bits/stdc++.h>

using namespace std;
using ll = long long;

const int N = 1e5 + 5; // 示例数组大小
ll a[N]; // 原始数组

struct Node {
    ll g; // 区间 GCD
} tree[N << 2]; // 线段树数组，大小开 4 倍

ll gcd(ll x, ll y) { // 保证处理 0 的情况
    if (x == 0) return y;
    if (y == 0) return x;
    while (y != 0) {
        ll r = x % y;
        x = y;
        y = r;
    }
    return x;
}

void pushup(int p) {
    tree[p].g = gcd(tree[p << 1].g, tree[p << 1 | 1].g);
}

void build(int p, int l, int r) {
    if (l == r) {
        tree[p].g = a[l]; // 叶子节点
        return;
    }
    int mid = (l + r) >> 1;
    build(p << 1, l, mid);
    build(p << 1 | 1, mid + 1, r);
    pushup(p);
}

// 查询区间 [ql, qr] 的 GCD
ll query(int p, int l, int r, int ql, int qr) {
    if (ql <= l && r <= qr) {
        return tree[p].g; // 完全覆盖
    }
    int mid = (l + r) >> 1;
    ll res_l = 0, res_r = 0; // 初始化为 0，因为 gcd(x, 0) = x
    if (ql <= mid) {
        res_l = query(p << 1, l, mid, ql, qr);
    }
    if (qr > mid) {
        res_r = query(p << 1 | 1, mid + 1, r, ql, qr);
    }
    return gcd(res_l, res_r); // 合并左右子区间结果
}

// 单点修改: 将 a[idx] 修改为 val
void update(int p, int l, int r, int idx, ll val) {
    if (l == r) {
        tree[p].g = val;
        a[idx] = val; // 更新原始数组 (如果需要)
        return;
    }
    int mid = (l + r) >> 1;
    if (idx <= mid) {
        update(p << 1, l, mid, idx, val);
    } else {
        update(p << 1 | 1, mid + 1, r, idx, val);
    }
    pushup(p);
}


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

    int n, q; // n: 数组大小, q: 查询次数
    cin >> n >> q;
    for (int i = 1; i <= n; ++i) {
        cin >> a[i];
    }

    build(1, 1, n); // 建树

    while (q--) {
        int type;
        cin >> type;
        if (type == 1) { // 查询操作
            int l, r;
            cin >> l >> r;
            cout << query(1, 1, n, l, r) << "\n";
        } else { // 修改操作 (示例)
            int idx;
            ll val;
            cin >> idx >> val;
            update(1, 1, n, idx, val);
        }
    }

    return 0;
}

稀疏表 (ST 表) 维护区间 GCD: ST 表适用于静态数组（无修改操作）的区间查询。预处理复杂度 O(n log n)，查询复杂度 O(1)。

• st[i][j] 表示区间 [i, i + 2^j - 1] 的 GCD。
• 预处理: st[i][0] = a[i]。st[i][j] = gcd(st[i][j-1], st[i + (1 << (j-1))][j-1])。
• 查询: 对于区间 [l, r]，长度 len = r - l + 1。令 k = log2(len)。查询结果为 gcd(st[l][k], st[r - (1 << k) + 1][k])。

ST 表在纯查询场景下效率更高。

6. 竞赛中的注意事项与技巧

• 数据范围与溢出: GCD 计算本身不易溢出，但涉及 LCM 或 exGCD 解乘法时，中间结果或最终解可能非常大，务必使用 long long。特别注意 LCM 的计算方式。
• 效率: O(log N) 的复杂度通常足够快。但在极端情况下（如在复杂度较高的算法内部频繁调用 GCD），需要考虑常数优化（迭代代替递归）或是否有更优的整体思路。
• GCD 的性质应用: 灵活运用 gcd(a, b) = gcd(a, b-a)、gcd(ka, kb) = k*gcd(a, b) 等性质有时能简化问题。例如，gcd(a, b) = gcd(a, a+b)。
• 与互质的关系: gcd(a, b) = d 等价于 a = d*a', b = d*b' 且 gcd(a', b') = 1。这个转换在很多数论证明和解题中非常有用。
• 编码细节:
  • 处理 gcd(a, 0) 的情况。标准欧几里得算法实现能正确处理。
  • exGCD 求解 ax + by = c 时，注意 c 必须是 gcd(a, b) 的倍数。
  • 求模逆元时，确保 gcd(a, m) = 1。得到的解 x 需要调整到 [0, m-1] 范围内，使用 (x % m + m) % m。

7. 总结

最大公约数 GCD 是算法竞赛数论部分的基础工具。掌握欧几里得算法（递归与迭代）、扩展欧几里得算法及其在线性丢番图方程、模逆元求解中的应用是基本要求。理解 GCD 与 LCM 的关系、多项 GCD/LCM 的计算、以及如何与数据结构结合处理区间 GCD 问题，能让你在面对更复杂的数论或组合问题时游刃有余。

8. 深入探索与进阶应用

在前文中，我们已经覆盖了 GCD 的核心算法和基本应用。现在，让我们挖掘一些更深的理论联系和更复杂的应用场景。

8.1 GCD 与素数分解 (Prime Factorization)

每个正整数 n > 1 都可以唯一地表示为其素因子的乘积（算术基本定理）：

$$n = p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}$$

其中 p_i 是不同的素数，a_i >= 1 是对应的指数。

GCD 和 LCM 可以通过素数分解来理解： 假设

$$a = p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}$$ $$b = p_1^{b_1} p_2^{b_2} \cdots p_k^{b_k}$$

（这里我们允许指数 a_i 或 b_i 为 0，以包含所有出现在 a 或 b 中的素因子 p_i）

那么：

$$gcd(a, b) = p_1^{\min(a_1, b_1)} p_2^{\min(a_2, b_2)} \cdots p_k^{\min(a_k, b_k)} $$$$lcm(a, b) = p_1^{\max(a_1, b_1)} p_2^{\max(a_2, b_2)} \cdots p_k^{\max(a_k, b_k)} $$

核心思想:

• GCD 取每个公共素因子的 最低 次幂。
• LCM 取每个出现过的素因子的 最高 次幂。

推论: 从这个角度可以轻松证明 gcd(a, b) * lcm(a, b) = a * b。 因为对于每个素因子 p_i，其在 gcd(a, b) 中的指数是 min(a_i, b_i)，在 lcm(a, b) 中的指数是 max(a_i, b_i)。 而在 a * b 中，p_i 的指数是 a_i + b_i。 由于 min(x, y) + max(x, y) = x + y，所以对于每个素因子 p_i，其在 gcd * lcm 中的指数等于在 a * b 中的指数。根据算术基本定理的唯一性，这两个数相等。

应用场景: 虽然直接用素数分解计算 GCD/LCM 通常比欧几里得算法慢（因为分解质因数困难），但这个理论视角在解决涉及多个数、需要分析因子结构或证明性质的问题时非常有用。例如，分析 gcd(a, lcm(b, c)) 或 lcm(a, gcd(b, c)) 等表达式。

8.2 二进制 GCD 算法 (Stein's Algorithm)

除了欧几里得算法，还有一种计算 GCD 的方法，特别是在不支持快速硬件除法或取模运算的机器上（或者在某些特定场景下，如处理大整数库内部实现时）可能更高效。它只使用减法、位移（除以2）和奇偶性判断。

算法原理基于以下观察:

1. 如果 a 和 b 都是偶数，gcd(a, b) = 2 * gcd(a/2, b/2)。
2. 如果 a 是偶数，b 是奇数，gcd(a, b) = gcd(a/2, b)。（因为 2 不可能是公约数）
3. 如果 a 是奇数，b 是偶数，gcd(a, b) = gcd(a, b/2)。
4. 如果 a 和 b 都是奇数，且 a > b，则 gcd(a, b) = gcd(a - b, b)。注意 a - b 是偶数。这步可以将问题转化为包含偶数的情况。
5. gcd(a, 0) = a。

实现:

#include <bits/stdc++.h>

using namespace std;
using ll = long long;

// 二进制 GCD 算法 (Stein's Algorithm)
// 时间复杂度: O(log(max(a, b)))， 位运算通常更快
ll binary_gcd(ll a, ll b) {
    if (a == 0) return b;
    if (b == 0) return a;
    // 处理负数
    a = abs(a);
    b = abs(b);

    // 计算 a 和 b 中因子 2 的总次数 k
    int k = 0;
    while (((a | b) & 1) == 0) { // a 和 b 都是偶数
        a >>= 1; // a /= 2
        b >>= 1; // b /= 2
        k++;     // 记录公共因子 2 的数量
    }

    // 现在 a 和 b 中至少有一个是奇数
    while ((a & 1) == 0) { // 如果 a 是偶数
        a >>= 1;
    }

    // 现在 a 肯定是奇数
    do {
        while ((b & 1) == 0) { // 如果 b 是偶数
            b >>= 1;
        }
        // 现在 a 和 b 都是奇数
        if (a > b) {
            swap(a, b); // 保证 a <= b
        }
        b = (b - a); // b 变为偶数，下次循环 b 会被除 2
    } while (b != 0);

    // 结果是 a 乘以之前提取的公共因子 2^k
    return a << k; // a * (2^k)
}

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

    ll a, b;
    cin >> a >> b;
    cout << binary_gcd(a, b) << endl;
    // cout << std::gcd(a, b) << endl; // 对比 STL 实现

    return 0;
}

复杂度与优势:

• 时间复杂度: 仍然是 O(log(max(a, b)))。每次操作（除以2或相减）都会使至少一个数减小。
• 优势: 避免了取模和除法运算，只用位运算和减法。在某些底层实现或特定硬件上可能更快。理解其原理也有助于拓宽思路。
• 劣势: 代码相对欧几里得算法略长。在有快速硬件除法/取模的现代 CPU 上，实际性能优势通常不明显，有时甚至可能因为分支预测等因素稍慢。

8.3 拉梅定理 (Lamé's Theorem) 与复杂度界限

欧几里得算法的效率非常高，但其精确的步数上界是什么？拉梅定理给出了答案。

定理: 对于两个正整数 a 和 b (a > b)，使用欧几里得算法计算 gcd(a, b) 所需的除法（或取模）次数不超过 b 的十进制位数的 5 倍。 更精确的界限与斐波那契数列有关：如果欧几里得算法需要 k 步，那么 a >= F_{k+2} 且 b >= F_{k+1}，其中 F_n 是第 n 个斐波那契数（F_1=1, F_2=1, ...）。 由于斐波那契数大致按黄金分割比 phi = (1+sqrt(5))/2 ≈ 1.618 的幂增长，即 F_n ≈ phi^n / sqrt(5)，这表明步数 k 大约是 log_phi(b)，从而证明了 O(log b) 的复杂度。

意义: 拉梅定理为欧几里得算法提供了一个坚实的理论基础，证明了它的效率，并揭示了其最坏情况与斐波那契数列的深刻联系（当输入是相邻的斐波那契数时，算法步数最多）。

8.4 GCD 与欧拉函数 (Euler's Totient Function)

欧拉函数 phi(n) 定义为小于或等于 n 的正整数中与 n 互质（即 GCD 为 1）的数的个数。 GCD 和 phi(n) 之间有许多深刻的联系。

关键性质:

1. 高斯定理 (Gauss's Theorem):

   $$\sum_{d|n} \phi(d) = n$$

   其中求和遍历 n 的所有正约数 d。这个性质非常重要，经常用于涉及 GCD 的求和问题。

   • 直观理解: 考虑分数 1/n, 2/n, ..., n/n。将它们化为最简分数 k/d，其中 d 必须是 n 的约数。对于每个约数 d，分母为 d 的最简分数恰好有 phi(d) 个（分子 k 满足 1 <= k <= d 且 gcd(k, d) = 1）。所有这些最简分数加起来正好是 n 个原始分数。

2. GCD 的分布: 对于给定的 n，考虑 gcd(1, n), gcd(2, n), ..., gcd(n, n)。这些 GCD 值有什么分布规律？ 一个数 k (1 <= k <= n) 满足 gcd(k, n) = d 的充要条件是：

   • d 必须是 n 的约数 (d | n)。
   • d | k。
   • gcd(k/d, n/d) = 1。

   令 k = d * i，则 1 <= d*i <= n，即 1 <= i <= n/d。 条件变为 gcd(i, n/d) = 1。 满足 1 <= i <= n/d 且 gcd(i, n/d) = 1 的 i 的个数正好是 phi(n/d)。 因此，对于 n 的任意约数 d，满足 gcd(k, n) = d 的 k (1 <= k <= n) 的个数为 phi(n/d)。 这也验证了高斯定理：sum_{d|n} phi(n/d) = sum_{d'|n} phi(d') = n （因为 d 遍历 n 的约数时，n/d 也遍历 n 的约数）。

3. 另一个恒等式:

   $$gcd(a, b) = \sum_{d|a, d|b} \phi(d)$$

   这个恒等式在理论推导中可能出现，但直接计算通常不实用。

应用：GCD 求和问题 一个经典问题是计算 sum_{k=1}^n gcd(k, n)。 利用上面的分布规律：

$$\sum_{k=1}^{n} gcd(k, n) = \sum_{d|n} \sum_{k=1, gcd(k,n)=d}^{n} d $$$$= \sum_{d|n} d \times (\text{number of } k \in [1, n] \text{ such that } gcd(k, n) = d) $$$$= \sum_{d|n} d \times \phi(n/d)$$

这个和可以通过枚举 n 的所有约数 d 来计算。如果 n 很大但其约数个数不多（例如，n 是素数的幂或少量素数的乘积），或者可以快速分解质因数，这种方法是可行的。计算 phi 值可以基于 n 的质因数分解。

代码示例：计算 sum_{k=1}^n gcd(k, n)

#include <bits/stdc++.h>

using namespace std;
using ll = long long;

// 计算 phi(n)
// 需要先分解质因数，或者直接用公式
// 如果需要多次计算 phi，可以预处理素数筛
ll phi(ll n) {
    ll res = n;
    for (ll i = 2; i * i <= n; ++i) {
        if (n % i == 0) {
            res = res / i * (i - 1);
            while (n % i == 0) {
                n /= i;
            }
        }
    }
    if (n > 1) { // 处理最后一个大的质因子
        res = res / n * (n - 1);
    }
    return res;
}

// 欧几里得算法 (如果下面直接用 gcd 计算对比)
ll gcd(ll a, ll b) {
    return b == 0 ? a : gcd(b, a % b);
}


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

    ll n;
    cin >> n;

    ll sum_gcd = 0;
    // 方法一: 利用公式 sum_{d|n} d * phi(n/d)
    for (ll d = 1; d * d <= n; ++d) {
        if (n % d == 0) {
            sum_gcd += d * phi(n / d);
            if (d * d != n) { // 避免重复计算完全平方数
                ll d2 = n / d;
                sum_gcd += d2 * phi(n / d2);
            }
        }
    }
    cout << "Sum gcd(k, n) [Formula]: " << sum_gcd << endl;

    // 方法二: 直接暴力计算 (用于小 n 验证)
    ll direct_sum = 0;
    for (ll k = 1; k <= n; ++k) {
        direct_sum += gcd(k, n);
    }
     cout << "Sum gcd(k, n) [Direct]: " << direct_sum << endl;


    return 0;
}

更复杂的问题: 类似地，可以推广到计算 sum_{i=1}^n sum_{j=1}^m f(gcd(i, j)) 这类二维 GCD 求和问题。这类问题通常需要更高级的技术，如莫比乌斯反演 (Möbius Inversion) 或基于 sum_{d|n} phi(d) = n 的变换技巧。

8.5 几何意义：格点可见性 (Lattice Point Visibility)

考虑二维平面上的整点（格点）(x, y)，其中 x, y 是整数。 一个格点 (x, y) 被称为从原点 (0, 0) 可见的，如果连接 (0, 0) 和 (x, y) 的线段上没有其他格点（除了端点）。

定理: 格点 (x, y) (不包括原点) 从原点可见，当且仅当 gcd(|x|, |y|) = 1。

证明思路:

• 如果 d = gcd(|x|, |y|) > 1，则 x = d*x'，y = d*y'。那么格点 (x', y') 就在 (0, 0) 和 (x, y) 之间的线段上，且 (x', y') 不是端点（因为 d > 1）。所以 (x, y) 不可见。
• 如果 gcd(|x|, |y|) = 1，假设线段上存在另一个格点 (x', y')。那么 (x', y') 必须是 (x, y) 的一个比例缩放，即 x' = kx, y' = ky 对于某个 0 < k < 1。由于 x', y' 是整数，k 必须是有理数，设 k = p/q（最简形式）。则 x' = px/q, y' = py/q。由于 x', y' 是整数，q 必须整除 px 和 py。因为 p, q 互质，所以 q 必须整除 x 和 y。但已知 gcd(x, y) = 1，所以 q 只能是 1。这与 k < 1 矛盾。因此线段中间没有其他格点，(x, y) 可见。

应用: 这个问题引出了计数问题，例如：在一个 N x M 的矩形区域内（1 <= x <= N, 1 <= y <= M），有多少个格点从原点可见？这等价于计算满足 1 <= x <= N, 1 <= y <= M 且 gcd(x, y) = 1 的数对 (x, y) 的数量。这类问题通常也与欧拉函数或莫比乌斯反演有关。

8.6 多项式 GCD (Polynomial GCD)

GCD 的概念可以推广到多项式环。对于两个系数在某个域（例如有理数、实数、复数或有限域）上的多项式 P(x) 和 Q(x)，它们的最大公约式 gcd(P(x), Q(x)) 是能同时整除 P(x) 和 Q(x) 的次数最高的多项式（通常要求首项系数为 1 或某个标准形式）。

计算多项式 GCD 同样可以使用欧几里得算法，只需将整数的除法和取模替换为多项式的带余除法。

应用:

• 求多项式的公共根。r 是 P(x) 和 Q(x) 的公共根当且仅当 (x-r) 是 gcd(P(x), Q(x)) 的因子。
• 化简有理函数 P(x) / Q(x)。可以约掉公因子 gcd(P(x), Q(x))。
• 在某些编码理论或代数问题中。