跳转至

# 最大公约数

## 定义

最大公约数即为 Greatest Common Divisor，常缩写为 gcd．

一组整数的公约数，是指同时是这组数中每一个数的约数的数．±1±1 是任意一组整数的公约数．

一组整数的最大公约数，是指所有公约数里面最大的一个．

对不全为 0 0 的整数 𝑎,𝑏 a,b，将其最大公约数记为 g c d(𝑎,𝑏)gcd(a,b)，不引起歧义时可简写为 (𝑎,𝑏)(a,b)．

对不全为 0 0 的整数 𝑎 1,…,𝑎 𝑛 a 1,…,a n，将其最大公约数记为 g c d(𝑎 1,…,𝑎 𝑛)gcd(a 1,…,a n)，不引起歧义时可简写为 (𝑎 1,…,𝑎 𝑛)(a 1,…,a n)．

最大公约数与最小公倍数的性质见 数论基础．

那么如何求最大公约数呢？我们先考虑两个数的情况．

### 欧几里得算法

#### 过程

如果我们已知两个数 𝑎 a 和 𝑏 b，如何求出二者的最大公约数呢？

不妨设 𝑎 >𝑏 a>b．

我们发现如果 𝑏 b 是 𝑎 a 的约数，那么 𝑏 b 就是二者的最大公约数． 下面讨论不能整除的情况，即 𝑎 =𝑏 ×𝑞 +𝑟 a=b×q+r，其中 𝑟 <𝑏 r<b．

我们通过证明可以得到 g c d(𝑎,𝑏) =g c d(𝑏,𝑎 m o d 𝑏)gcd(a,b)=gcd(b,a mod b)，过程如下：

证明

设 𝑎 =𝑏 𝑘 +𝑐 a=b k+c，显然有 𝑐 =𝑎 m o d 𝑏 c=a mod b．设 𝑑 ∣𝑎, 𝑑 ∣𝑏 d∣a, d∣b，则 𝑐 =𝑎 −𝑏 𝑘,𝑐 𝑑 =𝑎 𝑑 −𝑏 𝑑 𝑘 c=a−b k,c d=a d−b d k．

由右边的式子可知 𝑐 𝑑 c d 为整数，即 𝑑 ∣𝑐 d∣c，所以对于 𝑎,𝑏 a,b 的公约数，它也会是 𝑏,𝑎 m o d 𝑏 b,a mod b 的公约数．

反过来也需要证明：

设 𝑑 ∣𝑏, 𝑑 ∣(𝑎 m o d 𝑏)d∣b, d∣(a mod b)，我们还是可以像之前一样得到以下式子 𝑎 m o d 𝑏 𝑑 =𝑎 𝑑 −𝑏 𝑑 𝑘, 𝑎 m o d 𝑏 𝑑 +𝑏 𝑑 𝑘 =𝑎 𝑑 a mod b d=a d−b d k, a mod b d+b d k=a d．

因为左边式子显然为整数，所以 𝑎 𝑑 a d 也为整数，即 𝑑 ∣𝑎 d∣a，所以 𝑏,𝑎 m o d 𝑏 b,a mod b 的公约数也是 𝑎,𝑏 a,b 的公约数．

既然两式公约数都是相同的，那么最大公约数也会相同．

所以得到式子 g c d(𝑎,𝑏) =g c d(𝑏,𝑎 m o d 𝑏)gcd(a,b)=gcd(b,a mod b)

既然得到了 g c d(𝑎,𝑏) =g c d(𝑏,𝑟)gcd(a,b)=gcd(b,r)，这里两个数的大小是不会增大的，那么我们也就得到了关于两个数的最大公约数的一个递归求法．

#### 实现

C++Java Python

    1
    2
    3
    4
    5
    6
    7
    8
    // Version 1
    int gcd(int a, int b) {
      if (b == 0) return a;
      return gcd(b, a % b);
    }

    // Version 2
    int gcd(int a, int b) { return b == 0 ? a : gcd(b, a % b); }

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    // Version 1
    public int gcd(int a, int b) {
        if (b == 0) return a;
        return gcd(b, a % b);
    }

    // Version 2
    public int gcd(int a, int b) {
        return b == 0 ? a : gcd(b, a % b);
    }

    1
    2
    3
    4
    def gcd(a, b):
        if b == 0:
            return a
        return gcd(b, a % b)

递归至 `b == 0`（即上一步的 `a % b == 0`）的情况再返回值即可．

根据上述递归求法，我们也可以写出一个迭代求法：

C++Java Python

    1
    2
    3
    4
    5
    6
    7
    8
    int gcd(int a, int b) {
      while (b != 0) {
        int tmp = a;
        a = b;
        b = tmp % b;
      }
      return a;
    }

    1
    2
    3
    4
    5
    6
    7
    8
    public int gcd(int a, int b) {
        while(b != 0) {
            int tmp = a;
            a = b;
            b = tmp % b;
        }
        return a;
    }

    1
    2
    3
    4
    def gcd(a, b):
        while b != 0:
            a, b = b, a % b
        return a

上述算法都可被称作欧几里得算法（Euclidean algorithm）．

另外，对于 C++17，我们可以使用 `<numeric>` 头中的 `std::gcd` 与 `std::lcm` 来求最大公约数和最小公倍数．

注意

在部分编译器中，C++14 中可以用 `std::__gcd(a,b)` 函数来求最大公约数，但是其仅作为 `std::rotate` 的私有辅助函数．1使用该函数可能会导致预期之外的问题，故一般情况下不推荐使用．

如果两个数 𝑎 a 和 𝑏 b 满足 g c d(𝑎,𝑏) =1 gcd(a,b)=1，我们称 𝑎 a 和 𝑏 b 互质．

#### 性质

欧几里得算法的时间效率如何呢？下面我们证明，在输入为两个长为 𝑛 n 的二进制整数时，欧几里得算法的时间复杂度为 𝑂(𝑛)O(n)．（换句话说，在默认 𝑎,𝑏 a,b 同阶的情况下，时间复杂度为 𝑂(l o g⁡m a x(𝑎,𝑏))O(log⁡max(a,b))．）

证明

当我们求 g c d(𝑎,𝑏)gcd(a,b) 的时候，会遇到两种情况：

- 𝑎 <𝑏 a<b，这时候 g c d(𝑎,𝑏) =g c d(𝑏,𝑎)gcd(a,b)=gcd(b,a)；
- 𝑎 ≥𝑏 a≥b，这时候 g c d(𝑎,𝑏) =g c d(𝑏,𝑎 m o d 𝑏)gcd(a,b)=gcd(b,a mod b)，而对 𝑎 a 取模会让 𝑎 a 至少折半．这意味着这一过程最多发生 𝑂(l o g⁡𝑎) =𝑂(𝑛)O(log⁡a)=O(n) 次．

第一种情况发生后一定会发生第二种情况，因此第一种情况的发生次数一定 不多于 第二种情况的发生次数．

从而我们最多递归 𝑂(𝑛)O(n) 次就可以得出结果．

事实上，假如我们试着用欧几里得算法去求 斐波那契数列 相邻两项的最大公约数，会让该算法达到最坏复杂度．

### 更相减损术

大整数取模的时间复杂度较高，而加减法时间复杂度较低．针对大整数，我们可以用加减代替乘除求出最大公约数．

#### 过程

已知两数 𝑎 a 和 𝑏 b，求 g c d(𝑎,𝑏)gcd(a,b)．

不妨设 𝑎 ≥𝑏 a≥b，若 𝑎 =𝑏 a=b，则 g c d(𝑎,𝑏) =𝑎 =𝑏 gcd(a,b)=a=b． 否则，∀𝑑 ∣𝑎,𝑑 ∣𝑏∀d∣a,d∣b，可以证明 𝑑 ∣𝑎 −𝑏 d∣a−b．

因此，𝑎 a 和 𝑏 b 的 所有 公因数都是 𝑎 −𝑏 a−b 和 𝑏 b 的公因数，g c d(𝑎,𝑏) =g c d(𝑎 −𝑏,𝑏)gcd(a,b)=gcd(a−b,b)．

#### Stein 算法的优化

如果 𝑎 ≫𝑏 a≫b，更相减损术的 𝑂(𝑛)O(n) 复杂度将会达到最坏情况．

考虑一个优化，若 2 ∣𝑎,2 ∣𝑏 2∣a,2∣b，g c d(𝑎,𝑏) =2 g c d(𝑎 2,𝑏 2)gcd(a,b)=2 gcd(a 2,b 2)．

否则，若 2 ∣𝑎 2∣a（2 ∣𝑏 2∣b 同理），因为 2 ∣𝑏 2∣b 的情况已经讨论过了，所以 2 ∤𝑏 2∤b．因此 g c d(𝑎,𝑏) =g c d(𝑎 2,𝑏)gcd(a,b)=gcd(a 2,b)．

优化后的算法（即 Stein 算法）时间复杂度是 𝑂(l o g⁡𝑛)O(log⁡n)．

证明

若 2 ∣𝑎 2∣a 或 2 ∣𝑏 2∣b，每次递归至少会将 𝑎,𝑏 a,b 之一减半．

否则，2 ∣𝑎 −𝑏 2∣a−b，回到了上一种情况．

算法最多递归 𝑂(l o g⁡𝑛)O(log⁡n) 次．

#### 实现

高精度模板见 高精度计算．

高精度运算需实现：减法、大小比较、左移、右移（可用低精乘除代替）、二进制末位 0 的个数（可以通过判断奇偶暴力计算）．

C++

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    Big gcd(Big a, Big b) {
      if (a == 0) return b;
      if (b == 0) return a;
      // 记录a和b的公因数2出现次数，countr_zero表示二进制末位0的个数
      int atimes = countr_zero(a);
      int btimes = countr_zero(b);
      int mintimes = min(atimes, btimes);
      a >>= atimes;
      for (;;) {
        // a和b公因数中的2已经计算过了，后面不可能出现a为偶数的情况
        b >>= btimes;
        // 确保 a<=b
        if (a > b) swap(a, b);
        b -= a;
        if (b == 0) break;
        btimes = countr_zero(b);
      }
      return a << mintimes;
    }

上述代码参考了 libstdc++ 和 MSVC 对 C++17 `std::gcd` 的实现．在 `unsigned int` 和 `unsigned long long` 的数据范围下，如果可以以极快的速度计算 `countr_zero`，则 Stein 算法比欧几里得算法来得快，但反之则可能比欧几里得算法慢．

关于 countr_zero

1. gcc 有 内建函数 `__builtin_ctz`（32 位）或 `__builtin_ctzll`（64 位）可替换上述代码的 `countr_zero`；
2. 从 C++20 开始，头文件 `<bit>` 包含了 `std::countr_zero`；
3. 如果不使用不在标准库的函数，又无法使用 C++20 标准，下面的代码是一种在 Word-RAM with multiplication 模型下经过预处理后 𝑂(1)O(1) 的实现：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    constexpr int loghash[64] = {0,  32, 48, 56, 60, 62, 63, 31, 47, 55, 59, 61, 30,
                                 15, 39, 51, 57, 28, 46, 23, 43, 53, 58, 29, 14, 7,
                                 35, 49, 24, 44, 54, 27, 45, 22, 11, 37, 50, 25, 12,
                                 38, 19, 41, 52, 26, 13, 6,  3,  33, 16, 40, 20, 42,
                                 21, 10, 5,  34, 17, 8,  36, 18, 9,  4,  2,  1};

    int countr_zero(unsigned long long x) {
      return loghash[(x & -x) * 0x9150D32D8EB9EFC0Ui64 >> 58];
    }

而对于高精度运算，如果实现方法类似 `bitset`，则搭配上述对 `countr_zero` 的实现可以在 `O(n / w)` 的时间复杂度下完成．但如果不便按二进制位拆分，则只能暴力判断最大的 2 2 的幂因子，时间复杂度取决于实现．比如：

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    // 以小端序实现的二进制 Big，要求能枚举每一个元素
    int countr_zero(Big a) {
      int ans = 0;
      for (auto x : a) {
        if (x != 0) {
          ans += 32;  // 每一位数据类型的位长
        } else {
          return ans + countr_zero(x);
        }
      }
      return ans;
    }

    // 暴力计算，如需使用建议直接写进 gcd 加快常数
    int countr_zero(Big a) {
      int ans = 0;
      while ((a & 1) == 0) {
        a >>= 1;
        ++ans;
      }
      return ans;
    }

更多关于 `gcd` 实现上快慢的讨论可阅读 Fastest way to compute the greatest common divisor．

### 多个数的最大公约数

那怎么求多个数的最大公约数呢？显然答案一定是每个数的约数，那么也一定是每相邻两个数的约数．我们采用归纳法，可以证明，每次取出两个数求出答案后再放回去，不会对所需要的答案造成影响．

## 最小公倍数

接下来我们介绍如何求解最小公倍数（Least Common Multiple, LCM）．

### 定义

一组整数的公倍数，是指同时是这组数中每一个数的倍数的数．0 是任意一组整数的公倍数．

一组整数的最小公倍数，是指所有正的公倍数里面，最小的一个数．

对整数 𝑎,𝑏 a,b，将其最小公倍数记为 l c m⁡(𝑎,𝑏)lcm⁡(a,b)，不引起歧义时可简写为 [𝑎,𝑏][a,b]．

对整数 𝑎 1,…,𝑎 𝑛 a 1,…,a n，将其最小公倍数记为 l c m⁡(𝑎 1,…,𝑎 𝑛)lcm⁡(a 1,…,a n)，不引起歧义时可简写为 [𝑎 1,…,𝑎 𝑛][a 1,…,a n]．

### 两个数

设 𝑎 =𝑝 𝑘 𝑎 1 1 𝑝 𝑘 𝑎 2 2⋯𝑝 𝑘 𝑎 𝑠 𝑠 a=p 1 k a 1 p 2 k a 2⋯p s k a s，𝑏 =𝑝 𝑘 𝑏 1 1 𝑝 𝑘 𝑏 2 2⋯𝑝 𝑘 𝑏 𝑠 𝑠 b=p 1 k b 1 p 2 k b 2⋯p s k b s

我们发现，对于 𝑎 a 和 𝑏 b 的情况，二者的最大公约数等于

𝑝 m i n(𝑘 𝑎 1,𝑘 𝑏 1)1 𝑝 m i n(𝑘 𝑎 2,𝑘 𝑏 2)2⋯𝑝 m i n(𝑘 𝑎 𝑠,𝑘 𝑏 𝑠)𝑠 p 1 min(k a 1,k b 1)p 2 min(k a 2,k b 2)⋯p s min(k a s,k b s)

最小公倍数等于

𝑝 m a x(𝑘 𝑎 1,𝑘 𝑏 1)1 𝑝 m a x(𝑘 𝑎 2,𝑘 𝑏 2)2⋯𝑝 m a x(𝑘 𝑎 𝑠,𝑘 𝑏 𝑠)𝑠 p 1 max(k a 1,k b 1)p 2 max(k a 2,k b 2)⋯p s max(k a s,k b s)

由于 𝑘 𝑎 +𝑘 𝑏 =m a x(𝑘 𝑎,𝑘 𝑏) +m i n(𝑘 𝑎,𝑘 𝑏)k a+k b=max(k a,k b)+min(k a,k b)

所以得到结论是 g c d(𝑎,𝑏) ×l c m⁡(𝑎,𝑏) =𝑎 ×𝑏 gcd(a,b)×lcm⁡(a,b)=a×b

要求两个数的最小公倍数，先求出最大公约数即可．

### 多个数

可以发现，当我们求出两个数的 g c d gcd 时，求最小公倍数是 𝑂(1)O(1) 的复杂度．那么对于多个数，我们其实没有必要求一个共同的最大公约数再去处理，最直接的方法就是，当我们算出两个数的 g c d gcd，或许在求多个数的 g c d gcd 时候，我们将它放入序列对后面的数继续求解，那么，我们转换一下，直接将最小公倍数放入序列即可．

## 扩展欧几里得算法

扩展欧几里得算法（Extended Euclidean algorithm, EXGCD），常用于求 𝑎 𝑥 +𝑏 𝑦 =g c d(𝑎,𝑏)a x+b y=gcd(a,b) 的一组可行解．

### 过程

设

𝑎 𝑥 1 +𝑏 𝑦 1 =g c d(𝑎,𝑏)a x 1+b y 1=gcd(a,b)

𝑏 𝑥 2 +(𝑎 m o d 𝑏)𝑦 2 =g c d(𝑏,𝑎 m o d 𝑏)b x 2+(a mod b)y 2=gcd(b,a mod b)

由欧几里得定理可知：g c d(𝑎,𝑏) =g c d(𝑏,𝑎 m o d 𝑏)gcd(a,b)=gcd(b,a mod b)

所以 𝑎 𝑥 1 +𝑏 𝑦 1 =𝑏 𝑥 2 +(𝑎 m o d 𝑏)𝑦 2 a x 1+b y 1=b x 2+(a mod b)y 2

又因为 𝑎 m o d 𝑏 =𝑎 −(⌊𝑎 𝑏⌋ ×𝑏)a mod b=a−(⌊a b⌋×b)

所以 𝑎 𝑥 1 +𝑏 𝑦 1 =𝑏 𝑥 2 +(𝑎 −(⌊𝑎 𝑏⌋ ×𝑏))𝑦 2 a x 1+b y 1=b x 2+(a−(⌊a b⌋×b))y 2

𝑎 𝑥 1 +𝑏 𝑦 1 =𝑎 𝑦 2 +𝑏 𝑥 2 −⌊𝑎 𝑏⌋ ×𝑏 𝑦 2 =𝑎 𝑦 2 +𝑏(𝑥 2 −⌊𝑎 𝑏⌋𝑦 2)a x 1+b y 1=a y 2+b x 2−⌊a b⌋×b y 2=a y 2+b(x 2−⌊a b⌋y 2)

因为 𝑎 =𝑎,𝑏 =𝑏 a=a,b=b，所以 𝑥 1 =𝑦 2,𝑦 1 =𝑥 2 −⌊𝑎 𝑏⌋𝑦 2 x 1=y 2,y 1=x 2−⌊a b⌋y 2

将 𝑥 2,𝑦 2 x 2,y 2 不断代入递归求解直至 g c d gcd（最大公约数，下同）为 0 0 递归 𝑥 =1,𝑦 =0 x=1,y=0 回去求解．

### 实现

C++Python

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    int Exgcd(int a, int b, int &x, int &y) {
      if (!b) {
        x = 1;
        y = 0;
        return a;
      }
      int d = Exgcd(b, a % b, x, y);
      int t = x;
      x = y;
      y = t - (a / b) * y;
      return d;
    }

    1
    2
    3
    4
    5
    def Exgcd(a, b):
        if b == 0:
            return a, 1, 0
        d, x, y = Exgcd(b, a % b)
        return d, y, x - (a // b) * y

函数返回的值为 g c d gcd，在这个过程中计算 𝑥,𝑦 x,y 即可．

### 值域分析

𝑎 𝑥 +𝑏 𝑦 =g c d(𝑎,𝑏)a x+b y=gcd(a,b) 的解有无数个，显然其中有的解会爆 long long．
万幸的是，若 𝑏 ≠0 b≠0，扩展欧几里得算法求出的可行解必有 |𝑥| ≤𝑏,|𝑦| ≤𝑎|x|≤b,|y|≤a．
下面给出这一性质的证明．

证明

- g c d(𝑎,𝑏) =𝑏 gcd(a,b)=b 时，𝑎 m o d 𝑏 =0 a mod b=0，必在下一层终止递归．
  得到 𝑥 1 =0,𝑦 1 =1 x 1=0,y 1=1，显然 𝑎,𝑏 ≥1 ≥|𝑥 1|,|𝑦 1|a,b≥1≥|x 1|,|y 1|．
- g c d(𝑎,𝑏) ≠𝑏 gcd(a,b)≠b 时，设 |𝑥 2| ≤(𝑎 m o d 𝑏),|𝑦 2| ≤𝑏|x 2|≤(a mod b),|y 2|≤b．
  因为 𝑥 1 =𝑦 2,𝑦 1 =𝑥 2 −⌊𝑎 𝑏⌋𝑦 2 x 1=y 2,y 1=x 2−⌊a b⌋y 2
  所以 |𝑥 1| =|𝑦 2| ≤𝑏,|𝑦 1| ≤|𝑥 2| +|⌊𝑎 𝑏⌋𝑦 2| ≤(𝑎 m o d 𝑏) +⌊𝑎 𝑏⌋|𝑦 2||x 1|=|y 2|≤b,|y 1|≤|x 2|+|⌊a b⌋y 2|≤(a mod b)+⌊a b⌋|y 2|
  ≤𝑎 −⌊𝑎 𝑏⌋𝑏 +⌊𝑎 𝑏⌋|𝑦 2| ≤𝑎 −⌊𝑎 𝑏⌋(𝑏 −|𝑦 2|)≤a−⌊a b⌋b+⌊a b⌋|y 2|≤a−⌊a b⌋(b−|y 2|)
  𝑎 m o d 𝑏 =𝑎 −⌊𝑎 𝑏⌋𝑏 ≤𝑎 −⌊𝑎 𝑏⌋(𝑏 −|𝑦 2|) ≤𝑎 a mod b=a−⌊a b⌋b≤a−⌊a b⌋(b−|y 2|)≤a
  因此 |𝑥 1| ≤𝑏,|𝑦 1| ≤𝑎|x 1|≤b,|y 1|≤a 成立．

### 迭代法编写扩展欧几里得算法

首先，当 𝑥 =1 x=1，𝑦 =0 y=0，𝑥 1 =0 x 1=0，𝑦 1 =1 y 1=1 时，显然有：

{𝑎 𝑥+𝑏 𝑦=𝑎 𝑎 𝑥 1+𝑏 𝑦 1=𝑏{a x+b y=a a x 1+b y 1=b

成立．

已知 𝑎 m o d 𝑏 =𝑎 −(⌊𝑎 𝑏⌋ ×𝑏)a mod b=a−(⌊a b⌋×b)，下面令 𝑞 =⌊𝑎 𝑏⌋q=⌊a b⌋．参考迭代法求 gcd，每一轮的迭代过程可以表示为：

(𝑎,𝑏)→(𝑏,𝑎−𝑞 𝑏)(a,b)→(b,a−q b)

将迭代过程中的 𝑎 a 替换为 𝑎 𝑥 +𝑏 𝑦 =𝑎 a x+b y=a，𝑏 b 替换为 𝑎 𝑥 1 +𝑏 𝑦 1 =𝑏 a x 1+b y 1=b，可以得到：

{𝑎 𝑥+𝑏 𝑦=𝑎 𝑎 𝑥 1+𝑏 𝑦 1=𝑏→{𝑎 𝑥 1+𝑏 𝑦 1=𝑏 𝑎(𝑥−𝑞 𝑥 1)+𝑏(𝑦−𝑞 𝑦 1)=𝑎−𝑞 𝑏{a x+b y=a a x 1+b y 1=b→{a x 1+b y 1=b a(x−q x 1)+b(y−q y 1)=a−q b

据此就可以得到迭代法求 exgcd．

因为迭代的方法避免了递归，所以代码运行速度将比递归代码快一点．

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    int gcd(int a, int b, int& x, int& y) {
      x = 1, y = 0;
      int x1 = 0, y1 = 1, a1 = a, b1 = b;
      while (b1) {
        int q = a1 / b1;
        tie(x, x1) = make_tuple(x1, x - q * x1);
        tie(y, y1) = make_tuple(y1, y - q * y1);
        tie(a1, b1) = make_tuple(b1, a1 - q * b1);
      }
      return a1;
    }

如果你仔细观察 𝑎 1 a 1 和 𝑏 1 b 1，你会发现，他们在迭代版本的欧几里德算法中取值完全相同，并且以下公式无论何时（在 while 循环之前和每次迭代结束时）都是成立的：𝑥 ⋅𝑎 +𝑦 ⋅𝑏 =𝑎 1 x⋅a+y⋅b=a 1 和 𝑥 1 ⋅𝑎 +𝑦 1 ⋅𝑏 =𝑏 1 x 1⋅a+y 1⋅b=b 1．因此，该算法肯定能正确计算出 g c d gcd．

最后我们知道 𝑎 1 a 1 就是要求的 g c d gcd，有 𝑥 ⋅𝑎 +𝑦 ⋅𝑏 =𝑔 x⋅a+y⋅b=g．

#### 矩阵的解释

对于正整数 𝑎 a 和 𝑏 b 的一次辗转相除即 g c d(𝑎,𝑏) =g c d(𝑏,𝑎 m o d 𝑏)gcd(a,b)=gcd(b,a mod b) 使用矩阵表示如

[𝑏 𝑎 m o d 𝑏]=[0 1 1−⌊𝑎/𝑏⌋][𝑎 𝑏][b a mod b]=[0 1 1−⌊a/b⌋][a b]

其中向下取整符号 ⌊𝑐⌋⌊c⌋ 表示不大于 𝑐 c 的最大整数．我们定义变换 [𝑎 𝑏] ↦[0 1 1−⌊𝑎/𝑏⌋][𝑎 𝑏][a b]↦[0 1 1−⌊a/b⌋][a b]．

易发现欧几里得算法即不停应用该变换，有

[g c d(𝑎,𝑏)0]=(⋯[0 1 1−⌊𝑎/𝑏⌋][1 0 0 1])[𝑎 𝑏][gcd(a,b)0]=(⋯[0 1 1−⌊a/b⌋][1 0 0 1])[a b]

令

[𝑥 1 𝑥 2 𝑥 3 𝑥 4]=⋯[0 1 1−⌊𝑎/𝑏⌋][1 0 0 1][x 1 x 2 x 3 x 4]=⋯[0 1 1−⌊a/b⌋][1 0 0 1]

那么

[g c d(𝑎,𝑏)0]=[𝑥 1 𝑥 2 𝑥 3 𝑥 4][𝑎 𝑏][gcd(a,b)0]=[x 1 x 2 x 3 x 4][a b]

满足 𝑎 ⋅𝑥 1 +𝑏 ⋅𝑥 2 =g c d(𝑎,𝑏)a⋅x 1+b⋅x 2=gcd(a,b) 即扩展欧几里得算法，注意在最后乘了一个单位矩阵不会影响结果，提示我们可以在开始时维护一个 2 ×2 2×2 的单位矩阵编写更简洁的迭代方法如

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    int exgcd(int a, int b, int &x, int &y) {
      int x1 = 1, x2 = 0, x3 = 0, x4 = 1;
      while (b != 0) {
        int c = a / b;
        std::tie(x1, x2, x3, x4, a, b) =
            std::make_tuple(x3, x4, x1 - x3 * c, x2 - x4 * c, b, a - b * c);
      }
      x = x1, y = x2;
      return a;
    }

这种表述相较于递归更简单．

## 应用

- 10104 - Euclid Problem
- GYM - (J) once upon a time
- UVa - 12775 - Gift Dilemma

## 参考资料与链接

1. libstdc++: std Namespace Reference ↩

> 本页面最近更新：2026/1/27 12:26:08，更新历史
> 发现错误？想一起完善？ 在 GitHub 上编辑此页！
> 本页面贡献者：Ir1d, Tiphereth-A, Xeonacid, Enter-tainer, hsfzLZH1, c-forrest, iamtwz, ksyx, MegaOwIer, sshwy, StudyingFather, 383494, i-yyi, LuoshuiTianyi, mgt, untitledunrevised, Yanjun-Zhao, Backl1ght, buggg-hfc, FinParker, gi-b716, Great-designer, hhc0001, hly1204, hsiviter, huaruoji, Koishilll, Marcythm, Menci, NachtgeistW, ouuan, PwzXxm, Qubik65536, shawlleyw, tder6, TOMWT-qwq, VaneHsiung, warzone-oier, WillHouMoe
> 本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用