跳转至

最大公約數

定義

最大公約數即為 Greatest Common Divisor，常縮寫為 gcd。

一組整數的公約數，是指同時是這組數中每一個數的約數的數。\(\pm 1\) 是任意一組整數的公約數。

一組整數的最大公約數，是指所有公約數里面最大的一個。

那麼如何求最大公約數呢？我們先考慮兩個數的情況。

歐幾里得算法

過程

如果我們已知兩個數 \(a\) 和 \(b\)，如何求出二者的最大公約數呢？

不妨設 \(a > b\)

我們發現如果 \(b\) 是 \(a\) 的約數，那麼 \(b\) 就是二者的最大公約數。 下面討論不能整除的情況，即 \(a = b \times q + r\)，其中 \(r < b\)。

我們通過證明可以得到 \(\gcd(a,b)=\gcd(b,a \bmod b)\)，過程如下：

證明

設 \(a=bk+c\)，顯然有 \(c=a \bmod b\)。設 \(d \mid a,~d \mid b\)，則 \(c=a-bk, \frac{c}{d}=\frac{a}{d}-\frac{b}{d}k\)。

由右邊的式子可知 \(\frac{c}{d}\) 為整數，即 \(d \mid c\)，所以對於 \(a,b\) 的公約數，它也會是 \(b,a \bmod b\) 的公約數。

反過來也需要證明：

設 \(d \mid b,~d\mid (a \bmod b)\)，我們還是可以像之前一樣得到以下式子 \(\frac{a\bmod b}{d}=\frac{a}{d}-\frac{b}{d}k,~\frac{a\bmod b}{d}+\frac{b}{d}k=\frac{a}{d}\)。

因為左邊式子顯然為整數，所以 \(\frac{a}{d}\) 也為整數，即 \(d \mid a\)，所以 \(b,a\bmod b\) 的公約數也是 \(a,b\) 的公約數。

既然兩式公約數都是相同的，那麼最大公約數也會相同。

所以得到式子 \(\gcd(a,b)=\gcd(b,a\bmod b)\)

既然得到了 \(\gcd(a, b) = \gcd(b, r)\)，這裏兩個數的大小是不會增大的，那麼我們也就得到了關於兩個數的最大公約數的一個遞歸求法。

實現

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）的情況再返回值即可。

根據上述遞歸求法，我們也可以寫出一個迭代求法：

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\) 滿足 \(\gcd(a, b) = 1\)，我們稱 \(a\) 和 \(b\) 互質。

性質

歐幾里得算法的時間效率如何呢？下面我們證明，在輸入為兩個長為 \(n\) 的二進制整數時，歐幾里得算法的時間複雜度為 \(O(n)\)。（換句話説，在默認 \(a, b\) 同階的情況下，時間複雜度為 \(O(\log\max(a, b))\)。）

證明

當我們求 \(\gcd(a,b)\) 的時候，會遇到兩種情況：

  • \(a < b\)，這時候 \(\gcd(a,b)=\gcd(b,a)\)；
  • \(a \geq b\)，這時候 \(\gcd(a,b)=\gcd(b,a \bmod b)\)，而對 \(a\) 取模會讓 \(a\) 至少折半。這意味着這一過程最多發生 \(O(\log a) = O(n)\) 次。

第一種情況發生後一定會發生第二種情況，因此第一種情況的發生次數一定 不多於 第二種情況的發生次數。

從而我們最多遞歸 \(O(n)\) 次就可以得出結果。

事實上，假如我們試着用歐幾里得算法去求 斐波那契數列 相鄰兩項的最大公約數，會讓該算法達到最壞複雜度。

更相減損術

大整數取模的時間複雜度較高，而加減法時間複雜度較低。針對大整數，我們可以用加減代替乘除求出最大公約數。

過程

已知兩數 \(a\) 和 \(b\)，求 \(\gcd(a,b)\)。

不妨設 \(a \ge b\)，若 \(a = b\)，則 \(\gcd(a,b)=a=b\)。 否則，\(\forall d\mid a, d\mid b\)，可以證明 \(d\mid a-b\)。

因此，\(a\) 和 \(b\) 的 所有 公因數都是 \(a-b\) 和 \(b\) 的公因數，\(\gcd(a,b) = \gcd(a-b, b)\)。

Stein 算法的優化

如果 \(a\gg b\)，更相減損術的 \(O(n)\) 複雜度將會達到最壞情況。

考慮一個優化，若 \(2\mid a,2\mid b\)，\(\gcd(a,b) = 2\gcd\left(\dfrac a2, \dfrac b2\right)\)。

否則，若 \(2\mid a\)（\(2\mid b\) 同理），因為 \(2\mid b\) 的情況已經討論過了，所以 \(2 \nmid b\)。因此 \(\gcd(a,b)=\gcd\left(\dfrac a2,b\right)\)。

優化後的算法（即 Stein 算法）時間複雜度是 \(O(\log n)\)。

證明

若 \(2\mid a\) 或 \(2\mid b\)，每次遞歸至少會將 \(a,b\) 之一減半。

否則，\(2\mid a-b\)，回到了上一種情況。

算法最多遞歸 \(O(\log n)\) 次。

實現

高精度模板見 高精度計算。

高精度運算需實現：減法、大小比較、左移、右移（可用低精乘除代替）、判斷奇偶。

C++ 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
Big gcd(Big a, Big b) {
  // 記錄a和b的公因數2出現次數
  int atimes = 0, btimes = 0;
  while (a % 2 == 0) {
    a >>= 1;
    atimes++;
  }
  while (b % 2 == 0) {
    b >>= 1;
    btimes++;
  }
  for (;;) {
    // a和b公因數中的2已經計算過了，後面不可能出現a,b均為偶數的情況
    while (a % 2 == 0) {
      a >>= 1;
    }
    while (b % 2 == 0) {
      b >>= 1;
    }
    if (a == b) break;
    // 確保 a>=b
    if (a < b) swap(a, b);
    a -= b;
  }
  return a << min(atimes, btimes);
}

多個數的最大公約數

那怎麼求多個數的最大公約數呢？顯然答案一定是每個數的約數，那麼也一定是每相鄰兩個數的約數。我們採用歸納法，可以證明，每次取出兩個數求出答案後再放回去，不會對所需要的答案造成影響。

最小公倍數

接下來我們介紹如何求解最小公倍數（Least Common Multiple, LCM）。

定義

一組整數的公倍數，是指同時是這組數中每一個數的倍數的數。0 是任意一組整數的公倍數。

一組整數的最小公倍數，是指所有正的公倍數里面，最小的一個數。

兩個數

設 \(a = p_1^{k_{a_1}}p_2^{k_{a_2}} \cdots p_s^{k_{a_s}}\)，\(b = p_1^{k_{b_1}}p_2^{k_{b_2}} \cdots p_s^{k_{b_s}}\)

我們發現，對於 \(a\) 和 \(b\) 的情況，二者的最大公約數等於

\(p_1^{\min(k_{a_1}, k_{b_1})}p_2^{\min(k_{a_2}, k_{b_2})} \cdots p_s^{\min(k_{a_s}, k_{b_s})}\)

最小公倍數等於

\(p_1^{\max(k_{a_1}, k_{b_1})}p_2^{\max(k_{a_2}, k_{b_2})} \cdots p_s^{\max(k_{a_s}, k_{b_s})}\)

由於 \(k_a + k_b = \max(k_a, k_b) + \min(k_a, k_b)\)

所以得到結論是 \(\gcd(a, b) \times \operatorname{lcm}(a, b) = a \times b\)

要求兩個數的最小公倍數，先求出最大公約數即可。

多個數

可以發現，當我們求出兩個數的 \(\gcd\) 時，求最小公倍數是 \(O(1)\) 的複雜度。那麼對於多個數，我們其實沒有必要求一個共同的最大公約數再去處理，最直接的方法就是，當我們算出兩個數的 \(\gcd\)，或許在求多個數的 \(\gcd\) 時候，我們將它放入序列對後面的數繼續求解，那麼，我們轉換一下，直接將最小公倍數放入序列即可。

擴展歐幾里得算法

擴展歐幾里得算法（Extended Euclidean algorithm, EXGCD），常用於求 \(ax+by=\gcd(a,b)\) 的一組可行解。

過程

設

\(ax_1+by_1=\gcd(a,b)\)

\(bx_2+(a\bmod b)y_2=\gcd(b,a\bmod b)\)

由歐幾里得定理可知：\(\gcd(a,b)=\gcd(b,a\bmod b)\)

所以 \(ax_1+by_1=bx_2+(a\bmod b)y_2\)

又因為 \(a\bmod b=a-(\lfloor\frac{a}{b}\rfloor\times b)\)

所以 \(ax_1+by_1=bx_2+(a-(\lfloor\frac{a}{b}\rfloor\times b))y_2\)

\(ax_1+by_1=ay_2+bx_2-\lfloor\frac{a}{b}\rfloor\times by_2=ay_2+b(x_2-\lfloor\frac{a}{b}\rfloor y_2)\)

因為 \(a=a,b=b\)，所以 \(x_1=y_2,y_1=x_2-\lfloor\frac{a}{b}\rfloor y_2\)

將 \(x_2,y_2\) 不斷代入遞歸求解直至 \(\gcd\)（最大公約數，下同）為 0 遞歸 x=1,y=0 回去求解。

實現

 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

函數返回的值為 \(\gcd\)，在這個過程中計算 \(x,y\) 即可。

值域分析

\(ax+by=\gcd(a,b)\) 的解有無數個，顯然其中有的解會爆 long long。
萬幸的是，若 \(b\not= 0\)，擴展歐幾里得算法求出的可行解必有 \(|x|\le b,|y|\le a\)。
下面給出這一性質的證明。

證明
  • \(\gcd(a,b)=b\) 時，\(a\bmod b=0\)，必在下一層終止遞歸。
    得到 \(x_1=0,y_1=1\)，顯然 \(a,b\ge 1\ge |x_1|,|y_1|\)。
  • \(\gcd(a,b)\not= b\) 時，設 \(|x_2|\le (a\bmod b),|y_2|\le b\)。
    因為 \(x_1=y_2,y_1=x_2-{\left\lfloor\dfrac{a}{b}\right\rfloor}y_2\)
    所以 \(|x_1|=|y_2|\le b,|y_1|\le|x_2|+|{\left\lfloor\dfrac{a}{b}\right\rfloor}y_2|\le (a\bmod b)+{\left\lfloor\dfrac{a}{b}\right\rfloor}|y_2|\)
    \(\le a-{\left\lfloor\dfrac{a}{b}\right\rfloor}b+{\left\lfloor\dfrac{a}{b}\right\rfloor}|y_2|\le a-{\left\lfloor\dfrac{a}{b}\right\rfloor}(b-|y_2|)\)
    \(a\bmod b=a-{\left\lfloor\dfrac{a}{b}\right\rfloor}b\le a-{\left\lfloor\dfrac{a}{b}\right\rfloor}(b-|y_2|)\le a\)
    因此 \(|x_1|\le b,|y_1|\le a\) 成立。

迭代法編寫擴展歐幾里得算法

首先，當 \(x = 1\)，\(y = 0\)，\(x_1 = 0\)，\(y_1 = 1\) 時，顯然有：

\[ \begin{cases} ax + by & = a \\ ax_1 + by_1 & = b \end{cases} \]

成立。

已知 \(a\bmod b = a - (\lfloor \frac{a}{b} \rfloor \times b)\)，下面令 \(q = \lfloor \frac{a}{b} \rfloor\)。參考迭代法求 gcd，每一輪的迭代過程可以表示為：

\[ (a, b) \rightarrow (b, a - qb) \]

將迭代過程中的 \(a\) 替換為 \(ax + by = a\)，\(b\) 替換為 \(ax_1 + by_1 = b\)，可以得到：

\[ \begin{aligned} & \begin{cases} ax + by & = a \\ ax_1 + by_1 & = b \end{cases} \\ \rightarrow & \begin{cases} ax_1 + by_1 & = b \\ a(x - qx_1) + b(y - qy_1) & = a - qb \end{cases} \end{aligned} \]

據此就可以得到迭代法求 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;
}

如果你仔細觀察 \(a_1\) 和 \(b_1\)，你會發現，他們在迭代版本的歐幾里德算法中取值完全相同，並且以下公式無論何時（在 while 循環之前和每次迭代結束時）都是成立的：\(x \cdot a +y \cdot b =a_1\) 和 \(x_1 \cdot a +y_1 \cdot b= b_1\)。因此，該算法肯定能正確計算出 \(\gcd\)。

最後我們知道 \(a_1\) 就是要求的 \(\gcd\)，有 \(x \cdot a +y \cdot b = g\)。

矩陣的解釋

對於正整數 \(a\) 和 \(b\) 的一次輾轉相除即 \(\gcd(a,b)=\gcd(b,a\bmod b)\) 使用矩陣表示如

\[ \begin{bmatrix} b\\a\bmod b \end{bmatrix} = \begin{bmatrix} 0&1\\1&-\lfloor a/b\rfloor \end{bmatrix} \begin{bmatrix} a\\b \end{bmatrix} \]

其中向下取整符號 \(\lfloor c\rfloor\) 表示不大於 \(c\) 的最大整數。我們定義變換 \(\begin{bmatrix}a\\b\end{bmatrix}\mapsto \begin{bmatrix}0&1\\1&-\lfloor a/b\rfloor\end{bmatrix}\begin{bmatrix}a\\b\end{bmatrix}\)。

易發現歐幾里得算法即不停應用該變換，有

\[ \begin{bmatrix} \gcd(a,b)\\0 \end{bmatrix} = \left( \cdots \begin{bmatrix} 0&1\\1&-\lfloor a/b\rfloor \end{bmatrix} \begin{bmatrix} 1&0\\0&1 \end{bmatrix} \right) \begin{bmatrix} a\\b \end{bmatrix} \]

令

\[ \begin{bmatrix} x_1&x_2\\x_3&x_4 \end{bmatrix} = \cdots \begin{bmatrix} 0&1\\1&-\lfloor a/b\rfloor \end{bmatrix} \begin{bmatrix} 1&0\\0&1 \end{bmatrix} \]

那麼

\[ \begin{bmatrix} \gcd(a,b)\\0 \end{bmatrix} = \begin{bmatrix} x_1&x_2\\x_3&x_4 \end{bmatrix} \begin{bmatrix} a\\b \end{bmatrix} \]

滿足 \(a\cdot x_1+b\cdot x_2=\gcd(a,b)\) 即擴展歐幾里得算法，注意在最後乘了一個單位矩陣不會影響結果，提示我們可以在開始時維護一個 \(2\times 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;
}

這種表述相較於遞歸更簡單。

應用

參考資料與鏈接

  1. libstdc++: std Namespace Reference ↩

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用