Skip to content

Greatest Common Divisor

Definition

The greatest common divisor (GCD) is often abbreviated as gcd.

A common divisor of a set of integers is a number that divides every number in the set. \(\pm 1\) is a common divisor of any set of integers.

The greatest common divisor of a set of integers is the largest among all common divisors.

For integers \(a,b\) not both zero, their greatest common divisor is denoted as \(\gcd(a,b)\), which can be abbreviated as \((a,b)\) when no ambiguity arises.

For integers \(a_1,\dots,a_n\) not all zero, their greatest common divisor is denoted as \(\gcd(a_1,\dots,a_n)\), which can be abbreviated as \((a_1,\dots,a_n)\) when no ambiguity arises.

For properties of the greatest common divisor and the least common multiple, see Number Theory Basics.

So how do we compute the greatest common divisor? Let's first consider the case of two numbers.

Euclidean Algorithm

Process

Given two numbers \(a\) and \(b\), how do we find their greatest common divisor?

Assume without loss of generality that \(a > b\).

If \(b\) divides \(a\), then \(b\) is the greatest common divisor of the two numbers. Now consider the case where \(b\) does not divide \(a\), i.e., \(a = b \times q + r\), where \(r < b\).

We can prove that \(\gcd(a,b)=\gcd(b,a \bmod b)\), as follows:

Proof

Let \(a=bk+c\), where clearly \(c=a \bmod b\). Let \(d \mid a,~d \mid b\), then \(c=a-bk, \frac{c}{d}=\frac{a}{d}-\frac{b}{d}k\).

From the right-hand side, \(\frac{c}{d}\) is an integer, so \(d \mid c\). Thus, any common divisor of \(a\) and \(b\) is also a common divisor of \(b\) and \(a \bmod b\).

Conversely, we also need to prove:

Let \(d \mid b,~d\mid (a \bmod b)\). We can derive the following as before: \(\frac{a\bmod b}{d}=\frac{a}{d}-\frac{b}{d}k,~\frac{a\bmod b}{d}+\frac{b}{d}k=\frac{a}{d}\).

Since the left-hand side is clearly an integer, \(\frac{a}{d}\) is also an integer, i.e., \(d \mid a\). Thus, any common divisor of \(b\) and \(a\bmod b\) is also a common divisor of \(a\) and \(b\).

Since the sets of common divisors are identical, the greatest common divisor must also be the same.

Hence, we obtain \(\gcd(a,b)=\gcd(b,a\bmod b)\).

Given \(\gcd(a, b) = \gcd(b, r)\), where the sizes of the two numbers do not increase, we have a recursive method to compute the greatest common divisor of two numbers.

Implementation

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)

Recursion stops when b == 0 (i.e., when a % b == 0 in the previous step) and returns the value.

Based on the recursive method above, we can also write an iterative version:

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

The above algorithms are known as the Euclidean algorithm.

Additionally, for C++17, we can use std::gcd and std::lcm from the <numeric> header to compute the greatest common divisor and least common multiple.

Note

In some compilers, the function std::__gcd(a,b) is available in C++14 to compute the greatest common divisor, but it is a private helper function for std::rotate.1 Using this function may lead to unexpected issues, so it is generally not recommended.

If two numbers \(a\) and \(b\) satisfy \(\gcd(a, b) = 1\), we say \(a\) and \(b\) are coprime.

Properties

What is the time efficiency of the Euclidean algorithm? We prove that when the input consists of two \(n\)-bit binary integers, the time complexity of the Euclidean algorithm is \(O(n)\). (In other words, when \(a\) and \(b\) are of the same order, the time complexity is \(O(\log\max(a, b))\).)

Proof

When computing \(\gcd(a,b)\), we encounter two cases:

  • \(a < b\), in which case \(\gcd(a,b)=\gcd(b,a)\);
  • \(a \geq b\), in which case \(\gcd(a,b)=\gcd(b,a \bmod b)\), and taking modulo \(a\) reduces \(a\) by at least half. This means this process occurs at most \(O(\log a) = O(n)\) times.

The first case always leads to the second case, so the number of occurrences of the first case is no more than that of the second case.

Thus, we recurse at most \(O(n)\) times to obtain the result.

In fact, if we try to compute the greatest common divisor of two consecutive Fibonacci numbers using the Euclidean algorithm, it will reach the worst-case complexity.

Subtraction-Based Euclidean Algorithm

Modulo operations on large integers are time-consuming, while addition and subtraction are faster. For large integers, we can use addition and subtraction instead of multiplication and division to compute the greatest common divisor.

Process

Given two numbers \(a\) and \(b\), find \(\gcd(a,b)\).

Assume without loss of generality that \(a \ge b\). If \(a = b\), then \(\gcd(a,b)=a=b\). Otherwise, \(\forall d\mid a, d\mid b\), it can be shown that \(d\mid a-b\).

Thus, all common divisors of \(a\) and \(b\) are also common divisors of \(a-b\) and \(b\), so \(\gcd(a,b) = \gcd(a-b, b)\).

Optimization via Stein's Algorithm

If \(a\gg b\), the \(O(n)\) complexity of the subtraction-based method will reach its worst case.

Consider an optimization: if \(2\mid a,2\mid b\), then \(\gcd(a,b) = 2\gcd\left(\dfrac a2, \dfrac b2\right)\).

Otherwise, if \(2\mid a\) (similarly for \(2\mid b\)), since the case \(2\mid b\) has already been discussed, we have \(2 \nmid b\). Thus, \(\gcd(a,b)=\gcd\left(\dfrac a2,b\right)\).

The optimized algorithm (Stein's algorithm) has a time complexity of \(O(\log n)\).

Proof

If \(2\mid a\) or \(2\mid b\), each recursion reduces at least one of \(a\) or \(b\) by half.

Otherwise, \(2\mid a-b\), which falls back to the previous case.

The algorithm recurses at most \(O(\log n)\) times.

Implementation

For high-precision templates, see High-Precision Computation.

High-precision operations require: subtraction, comparison, left shift, right shift (can be replaced by low-precision multiplication/division), and counting trailing zeros in binary (can be computed by checking parity).

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;
  // Record the number of times 2 appears as a common factor of a and b
  int atimes = countr_zero(a);
  int btimes = countr_zero(b);
  int mintimes = min(atimes, btimes);
  a >>= atimes;
  for (;;) {
    // The factor 2 in the common divisors of a and b has been accounted for; a cannot be even afterward.
    b >>= btimes;
    // Ensure a <= b
    if (a > b) swap(a, b);
    b -= a;
    if (b == 0) break;
    btimes = countr_zero(b);
  }
  return a << mintimes;
}

The above code references the implementations of std::gcd in libstdc++ and MSVC. For unsigned int and unsigned long long, if countr_zero can be computed extremely fast, Stein's algorithm is faster than the Euclidean algorithm; otherwise, it may be slower.

About countr_zero
  1. GCC has built-in functions __builtin_ctz (32-bit) or __builtin_ctzll (64-bit) that can replace countr_zero in the above code.
  2. Starting from C++20, the header <bit> includes std::countr_zero.
  3. If non-standard library functions are not allowed and C++20 is unavailable, the following code provides an \(O(1)\) implementation after preprocessing in the Word-RAM with multiplication model:
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];
}

For high-precision arithmetic, if the implementation is similar to bitset, combining the above countr_zero implementation can achieve \(O(n / w)\) time complexity. However, if binary splitting is not feasible, one must brute-force the largest power-of-two factor, with time complexity depending on the implementation. For example:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
// Binary Big implemented in little-endian, requiring enumeration of each element
int countr_zero(Big a) {
  int ans = 0;
  for (auto x : a) {
    if (x != 0) {
      ans += 32;  // Bit length of each data type
    } else {
      return ans + countr_zero(x);
    }
  }
  return ans;
}

// Brute-force computation; if used, it is recommended to embed directly into gcd for better performance.
int countr_zero(Big a) {
  int ans = 0;
  while ((a & 1) == 0) {
    a >>= 1;
    ++ans;
  }
  return ans;
}

For more discussion on the speed of gcd implementations, see Fastest way to compute the greatest common divisor.

Greatest Common Divisor of Multiple Numbers

How do we compute the greatest common divisor of multiple numbers? Clearly, the answer must be a divisor of each number and thus a divisor of every pair of adjacent numbers. By induction, we can prove that taking two numbers, computing their GCD, and putting the result back does not affect the final answer.

Least Common Multiple

Next, we introduce how to compute the least common multiple (LCM).

Definition

A common multiple of a set of integers is a number that is a multiple of every number in the set. 0 is a common multiple of any set of integers.

The least common multiple of a set of integers is the smallest positive common multiple.

For integers \(a,b\), their least common multiple is denoted as \(\operatorname{lcm}(a,b)\), which can be abbreviated as \([a,b]\) when no ambiguity arises.

For integers \(a_1,\dots,a_n\), their least common multiple is denoted as \(\operatorname{lcm}(a_1,\dots,a_n)\), which can be abbreviated as \([a_1,\dots,a_n]\) when no ambiguity arises.

Two Numbers

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

We observe that for \(a\) and \(b\), their greatest common divisor equals

\(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})}\)

and their least common multiple equals

\(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})}\)

Since \(k_a + k_b = \max(k_a, k_b) + \min(k_a, k_b)\),

we conclude that \(\gcd(a, b) \times \operatorname{lcm}(a, b) = a \times b\).

To compute the least common multiple of two numbers, first compute their greatest common divisor.

Multiple Numbers

We can see that once we compute the \(\gcd\) of two numbers, computing their LCM is \(O(1)\). For multiple numbers, there is no need to compute a common GCD first. Instead, when we compute the GCD of two numbers (or while computing the GCD of multiple numbers), we can replace them with their LCM in the sequence and continue.

Extended Euclidean Algorithm

The extended Euclidean algorithm (EXGCD) is commonly used to find a feasible solution to \(ax+by=\gcd(a,b)\).

Process

Let

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

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

By the Euclidean theorem: \(\gcd(a,b)=\gcd(b,a\bmod b)\)

Thus, \(ax_1+by_1=bx_2+(a\bmod b)y_2\)

Since \(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)\)

Since \(a=a,b=b\), we have \(x_1=y_2,y_1=x_2-\lfloor\frac{a}{b}\rfloor y_2\)

Recursively substitute \(x_2,y_2\) until \(\gcd\) (greatest common divisor) becomes \(0\), then backtrack with \(x=1,y=0\) to solve.

Implementation

 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

The function returns \(\gcd\), and during the process, \(x\) and \(y\) are computed.

Range Analysis

The equation \(ax+by=\gcd(a,b)\) has infinitely many solutions, and some may overflow long long.
Fortunately, if \(b\not= 0\), the feasible solution obtained by the extended Euclidean algorithm satisfies \(|x|\le b,|y|\le a\).
Below is a proof of this property.

Proof
  • When \(\gcd(a,b)=b\), \(a\bmod b=0\), and recursion terminates at the next level.
    We obtain \(x_1=0,y_1=1\), and clearly \(a,b\ge 1\ge |x_1|,|y_1|\).
  • When \(\gcd(a,b)\not= b\), assume \(|x_2|\le (a\bmod b),|y_2|\le b\).
    Since \(x_1=y_2,y_1=x_2-{\left\lfloor\dfrac{a}{b}\right\rfloor}y_2\),
    we have \(|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\)
    Thus, \(|x_1|\le b,|y_1|\le a\) holds.

Iterative Implementation of the Extended Euclidean Algorithm

First, when \(x = 1\), \(y = 0\), \(x_1 = 0\), \(y_1 = 1\), the following clearly holds:

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

Given \(a\bmod b = a - (\lfloor \frac{a}{b} \rfloor \times b)\), let \(q = \lfloor \frac{a}{b} \rfloor\). Referring to the iterative method for computing gcd, each iteration can be expressed as:

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

Replacing \(a\) with \(ax + by = a\) and \(b\) with \(ax_1 + by_1 = b\) in the iterative process, we get:

\[ \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} \]

This leads to the iterative method for exgcd.

Since iteration avoids recursion, the code runs slightly faster than the recursive version.

 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;
}

If you observe \(a_1\) and \(b_1\), you will notice they take exactly the same values as in the iterative Euclidean algorithm, and the following formulas hold at all times (before the while loop and at the end of each iteration): \(x \cdot a +y \cdot b =a_1\) and \(x_1 \cdot a +y_1 \cdot b= b_1\). Therefore, the algorithm correctly computes \(\gcd\).

Finally, we know \(a_1\) is the desired \(\gcd\), satisfying \(x \cdot a +y \cdot b = g\).

Matrix Interpretation

For one step of the Euclidean algorithm on positive integers \(a\) and \(b\), i.e., \(\gcd(a,b)=\gcd(b,a\bmod b)\), the matrix representation is:

\[ \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} \]

where the floor function \(\lfloor c\rfloor\) denotes the greatest integer less than or equal to \(c\). We define the transformation \(\begin{bmatrix}a\\b\end{bmatrix}\mapsto \begin{bmatrix}0&1\\1&-\lfloor a/b\rfloor\end{bmatrix}\begin{bmatrix}a\\b\end{bmatrix}\).

The Euclidean algorithm repeatedly applies this transformation:

\[ \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} \]

Let

\[ \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} \]

Then

\[ \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} \]

Satisfying \(a\cdot x_1+b\cdot x_2=\gcd(a,b)\) gives the extended Euclidean algorithm. Note that multiplying by the identity matrix at the end does not affect the result, suggesting we can maintain a \(2\times 2\) identity matrix initially for a more concise iterative implementation:

 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;
}

This representation is simpler than the recursive version.

Applications

  1. libstdc++: std Namespace Reference ↩

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