浅谈最大公约数

· · 算法·理论

写在前面

为备忘记录下，若要转载请标明出处（话说这么烂也没人会转载吧），感谢管理审核。

最大公约数

最大公约数（Greatest Common Divisor），简称 gcd。任意一组整数 a_1, a_2, a_3, ..., a_n 的公约数 d，当且仅当满足对于数列中的每一个数 a_i 有 d \mid a_i。而最大公约数则是所有可行的 d 中的最大的一个，记作 d = \gcd(a_1, a_2,..., a_n)。

性质

  1. 非负性，即对于任意整数 a_1, a_2, a_3, ..., a_{n - 1}, a_n，有 \gcd(a_1, a_2, ..., a_n) = \gcd(|a_1|, |a_2|, ..., |a_n|) \ge 0（由定义易证）。
  2. 轮换对称性（“交换性”似乎不太严谨），即 \gcd(a, b) = \gcd(b, a)（可以拓展到更多个数的情况，易证）。
  3. 结合性，即 \gcd(a, b, c) = \gcd(\gcd(a, b), c)（同样可以拓展到更多情况，易证）。
  5. 对于一组整数 a_1, a_2, ..., a_n，令 d = \gcd(a_1, a_2, ..., a_n)，则有 \gcd(\frac{a_1}{d}, \frac{a_2}{d}, ..., \frac{a_n}{d}) = 1，由定义易证。
  6. 随便根据定义推一下就有一大堆，再此不多赘述......

求法

这是本篇文章主要讨论的部分，下面会分别讲解一些常用的最大公约数求法。

欧几里得算法

欧几里得算法（Euclidean algorithm），又名辗转相除法，常用于求两数 a 和 b 的最大公约数。

算法原理

\gcd(a, b) = \gcd(b, a \bmod b)

证明

令 d = \gcd(a, b)，则有 a = x \times d， b = y \times d，x \nmid y。

令 k = \lfloor \frac{a}{b} \rfloor，那么有 a \bmod b = a - k \times b = x \times d - k \times y \times d = d \times (x - k \times y)。

又因为 x \nmid y，所以 \gcd(b, a \bmod b) = \gcd(y \times b, (x - k \times y) \times d) = d。

故 \gcd(a, b) = \gcd(b, a \bmod b)。

所以，我们可以递归地求得 \gcd(a, b)。当 b = 0 时，有上一层递归中的 a' \bmod b' = 0，即 b' = d，所以此时的 a 即为最大公约数，可以结束递归。

时间复杂度分析

若 a = b，一次递归后 b 变为 0，时间复杂度为 O(1)。

若 a > b，有如下情况:

若 a < b，一次递归后变为 \gcd(b, a \bmod b) = \gcd(b, a)，情况等同于 a' > b'，时间复杂度为 O(1 + \log \max(a, b)) \approx O(\log \max(a, b))。

所以，欧几里得算法的平均时间复杂度为 O(\log \max(a, b))，最坏时间复杂度为 O(\log \max(a, b))，一般来说是可以视为 O(1) 的。

Code

#include <bits/stdc++.h>
using namespace std;
int x, y;
int gcd(int a, int b){
    if (!b)return a;//若当前 b 等于 0，结束递归
    return gcd(b, a % b);
}
int main(){
    cin >> x >> y;
    cout << gcd(x, y);
    return 0;
}

更相减损术（Stein 优化）

欧几里得算法足以满足大部分的求解需求，但是取模运算的耗时较长，若在大整数运算中变会捉襟见肘。所以，一种新的解法应运而生，这便是更相减损术（Stein 优化）。

其运算过程为先判断参与运算的二者 a 和 b 的奇偶性，若都为偶数则先分出一个 2 并递归求解 \gcd(\frac{a}{2}, \frac{b}{2})，否则若只有一个偶数则递归求解 \gcd(\frac{a}{2}, b)，若全为奇数则递归求解 \gcd(b, \frac{a - b}{2})。

算法原理

  1. \gcd(a, b) = \gcd(b, a - b) \ (a \ge b)
  2. \gcd(a, b) = \gcd(\frac{a}{2}, b) \ (2 \mid a \text{ 且 } 2 \nmid b)

证明

原理 1：

令 \gcd(a, b) = d，则有 a = d \times x，b = d \times y，x \nmid y。

所以 a - b = d \times (x - y)，即 a 和 b 的所有公因数都是 b 和 a - b 的公因数。

继续证明 y \nmid x - y，令 b 和 a - b 的最大公因数为 h，则有 h \ge d 和 d \mid h，令 g = \frac{h}{d}，则有 b = n \times h = y \times g \times d 和 a - b = m \times h = (x - y) \times g \times d。

所以有 a = a - b + b = x \times g \times d = x \times d，即 g = 1，所以 h = d，可得 x \nmid x - y，得证。

原理 2（其实没什么证的）：

因为已知 2 \mid a 和 2 \nmid b，所以 2 \nmid \gcd(a, b)，那么就可以将 a 中的 2 分离出来，即 \gcd(a, b) = \gcd(\frac{a}{2}, b)，得证。

时间复杂度分析

综上，此算法的时间复杂度是是在 O(\log \max(a, b)) 量级的。

Code

#include <bits/stdc++.h>
using namespace std;
int x, y;
int gcd(int a, int b){
    if (!b)return a;
    if (b > a)swap(a, b);
    if (!(a & 1) && !(b & 1))return gcd(a >> 1, b >> 1) << 1;//都为偶数
    if ((a & 1) ^ (b & 1))return gcd(a & 1 ? a : a >> 1, b & 1 ? b : b >> 1);//一奇一偶
    return gcd(b, a - b >> 1);
}
int main(){
    cin >> x >> y;
    cout << gcd(x, y);
    return 0;
}

非递归版：

int gcd(int a, int b){
    if (a == b)return a;
    int az = ctz(a), az = ctz(b), z = az < bz ? az : bz, d;
    a >>= az;
    while (b)
        b >>= bz, bz = ctz(d = b - a),
        a = b < a ? b : a, b = abs(d);
    return a << z;
}//其中 ctz 为 __builtin_ctz，即数二进制末尾 0 的个数

拓展欧几里得算法

拓展欧几里得算法（Extended Euclidean algorithm, EXGCD），常被用于求解 ax + by = \gcd(a, b) 的一组可行解。

裴蜀定理

内容

若有整数 a、b，则一定存在整数 x、y 使得 ax + by = \gcd(a, b)。

即对于方程 ax + by = d \times \gcd(a, b)，若 a、b 为整数，则此方程有整数解。

证明

令集合 S = \{ax + by | x, y \in Z\}，记其中的最小正数为 m = ax_0 + by_0。

则对于任意 c \in S，有 c = qm + r \ (0 \le r < m \ | \ q \in Z) 和 c = ax_1+ by_1。

那么 r = c - qm = ax_1 + by_1 - aqx_0 - bqy_0 = a(x_1 - qx_0) + b(y_1 - qy_0) \in S。又因为 m 为 S 中最小正整数且 r \ge 0，所以 r = 0，即 m \mid c。

特别地，当 c 取 a 和 b 时，可得 m 为 a 和 b 的公约数。

设 n 为 a 和 b 的任意公约数，则有 n \mid m，也即 n \le m，所以 m 为 a 和 b 的最大公约数 \gcd(a, b)，得证。

实现

那么如何求 x 和 y 的具体值呢？这就是拓欧的天下了。

我们知道最大公约数可以使用辗转相除法求得，所以我们一步步由递归函数往上推。先把代码放上来：

int gcd(int a, int b){
    if (!b)return a;
    return gcd(b, a % b);
}

首先当 b' = 0 时，当前的 a' 即为 d = \gcd(a, b)，此时的 x' = 1, y' = 0。

那么 b' \ne 0 时，递归时则有

b'x'' + (a' \bmod b')y'' = \gcd(a', b')

也即

b'x'' + (a' - \left \lfloor \frac{a'}{b'} \right \rfloor b')y'' = \gcd(a', b')

化简得

a'y'' + b'(x'' - \left \lfloor \frac{a'}{b'} \right \rfloor y'') = \gcd(a', b')

所以当前层的 x' = y'', y' = x'' - \left \lfloor \frac{a'}{b'} \right \rfloor y''，即可递归求得最终的 x 和 y。

so the code：

void exgcd(int a, int b){
    if (!b){
        x = 1, y = 0;
        return;
    }
    exgcd(b, a % b);
    int now = x;
    x = y, y = now - a / b * y;
}

求通解

我们把 ax + by = \gcd(a, b) 拓展成 ax + by = c 的形式，考虑对其求解。

由裴蜀定理得若 \gcd(a, b) \nmid c 则方程无整数解，此处仅讨论有解的情况。

从 \gcd(a, b) 到 c 拓展了 \frac{c}{\gcd(a, b)} 倍，所以解也会对应地乘上 \frac{c}{\gcd(a, b)}。

即 x_0 = \frac{cx}{\gcd(a, b)}，y_0 = \frac{cy}{\gcd(a, b)}。

接下来求其通解，我们先让 x_1 = x_0 + n（n \in Z），看一下 y1 的变化。

由 ax_0 + by_0 = c 得 y_0 = \frac{c - ax_0}{b}，又有 ax_1 + by_1 = c 和 x_1 = x_0 + n，所以 y_1 = \frac{c - ax_1}{b} = \frac{c - a(x_0 + n)}{b} = \frac{c - ax_0}{b} - \frac{an}{b} = y_0 - \frac{a}{b}n。

把 \frac{a}{b}n 转化成可得到的最小整数，即同乘 \frac{b}{\gcd(a, b)} 得

x_1 = x_0 + \frac{b}{\gcd(a, b)}n y_1 = y_0 + \frac{a}{\gcd(a, b)}n

由此得到通解。

求最小非负整数解

那如何求最小非负整数解呢？

因为 x_1 = x_0 + \frac{b}{\gcd(a, b)}，所以必然存在使 x \equiv x_1 \pmod{\frac{b}{\gcd(a, b)}} 成立的解，所以可以得到 \min x_1 = \min (x_0 + \frac{b}{\gcd(a, b)}n) \bmod \frac{b}{\gcd(a, b)} = [(x_0 \bmod \frac{b}{\gcd(a, b)}) + \frac{b}{\gcd(a, b)}] \bmod \frac{b}{\gcd(a, b)}（x_1 \ge 0，加上 \frac{b}{\gcd(a, b)}) 是为防止负数）。

即

x_{min} = [(x_0 \bmod \frac{b}{\gcd(a, b)}) + \frac{b}{\gcd(a, b)}] \bmod \frac{b}{\gcd(a, b)}

同理得

y_{min} = [(y_0 \bmod \frac{a}{\gcd(a, b)}) + \frac{a}{\gcd(a, b)}] \bmod \frac{a}{\gcd(a, b)}

结语

乱写，有误还请指出，再次感谢审核。