AIOJ

Definition

The greatest common divisor is also called the Greatest Common Divisor and is commonly abbreviated as gcd.

A common divisor of a set of integers is a number that divides every integer in the set. $\pm 1$ is a common divisor of every set of integers.

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

For integers $0$ that are not both $a,b$ , their greatest common divisor is denoted by $\gcd(a,b)$ , or simply $(a,b)$ when there is no ambiguity.

For integers $0$ that are not all $a_1,\dots,a_n$ , their greatest common divisor is denoted by $\gcd(a_1,\dots,a_n)$ , or simply $(a_1,\dots,a_n)$ when there is no ambiguity.

See number theory basics for properties of the greatest common divisor and least common multiple.

How do we find a greatest common divisor? We first consider two numbers.

Euclidean algorithm

Process

Given two numbers $a$ and $b$ , how can we find their greatest common divisor?

Without loss of generality, assume $a > b$ .

If $b$ $a$ $b$ $a = b \times q + r$ $r < b$ divides , then is their greatest common divisor. Now consider the case where does not divide , namely with .

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

Proof

Let $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)$ , so clearly . If and , then and .

The right-hand side shows that is an integer, so . Therefore every common divisor of is also a common divisor of .

The converse must also be proved.

Let and . As above, we obtain and .

The left-hand side is clearly an integer, so is also an integer, hence . Thus every common divisor of is also a common divisor of .

Since the two pairs have exactly the same common divisors, their greatest common divisors are equal.

Therefore

Since $\gcd(a, b) = \gcd(b, r)$ and the two numbers do not increase, we obtain a recursive method for finding the greatest common divisor of two numbers.

Implementation

C++

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

Java

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

Python

def gcd(a, b):
    if b == 0:
        return a
    return gcd(b, a % b)

Return the value when the recursion reaches b == 0 (that is, when a % b == 0 in the previous step).

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

C++

int gcd(int a, int b) {
  while (b != 0) {
    int tmp = a;
    a = b;
    b = tmp % b;
  }
  return a;
}

Java

public int gcd(int a, int b) {
    while(b != 0) {
        int tmp = a;
        a = b;
        b = tmp % b;
    }
    return a;
}

Python

def gcd(a, b):
    while b != 0:
        a, b = b, a % b
    return a

All of the algorithms above are called the Euclidean algorithm.

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

Note

On some compilers, C++14 provides std::__gcd(a,b) for computing the greatest common divisor, but it is only a private helper function for std::rotate.¹ Using it may cause unexpected problems, so it is generally not recommended.

If two numbers $a$ $b$ $\gcd(a, b) = 1$ $a$ $b$ and satisfy , they are said to be coprime.

Properties

How efficient is the Euclidean algorithm? We prove that when the input consists of two $n$ $O(n)$ $a, b$ $O(\log\max(a, b))$ -bit binary integers, its time complexity is . (Equivalently, when and are assumed to have the same order of magnitude, the time complexity is .)

Proof

When computing $\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)$ , there are two cases:

, in which case ;
, in which case , and taking modulo at least halves . Thus this case occurs at most times.

The second case must occur after the first, so the first case occurs no more often than the second.

Consequently, at most recursive calls are needed.

In fact, applying the Euclidean algorithm to the greatest common divisor of two consecutive terms of the Fibonacci sequence makes the algorithm attain its worst-case complexity.

Subtractive Euclidean algorithm

Modulo operations on large integers have relatively high complexity, whereas addition and subtraction are cheaper. For large integers, we can use subtraction instead of multiplication and division to find the greatest common divisor.

Process

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

Assume $a \ge b$ $a = b$ $\gcd(a,b)=a=b$ $\forall d\mid a, d\mid b$ $d\mid a-b$ . If , then . Otherwise, for every such that and , we can prove that .

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

Stein algorithm optimization

If $a\gg b$ , the subtractive algorithm reaches its worst case with its $O(n)$ complexity.

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

Otherwise, suppose $2\mid a$ (the case $2\mid b$ is analogous). Since the case $2\mid b$ has already been considered, $2 \nmid b$ . Therefore $\gcd(a,b)=\gcd\left(\dfrac a2,b\right)$ .

The optimized algorithm, known as Stein's algorithm, has time complexity $O(\log n)$ .

Proof

If $2\mid a$ or $2\mid b$ , each recursive call at least halves one of $a,b$ .

Otherwise, $2\mid a-b$ , reducing the problem to the previous case.

The algorithm makes at most $O(\log n)$ recursive calls.

Implementation

See big integer arithmetic for a big-integer template.

Big-integer arithmetic must implement subtraction, comparison, left shift, right shift (which may be replaced by multiplication and division by small integers), and the number of trailing zero bits (which can be computed by checking parity directly).

C++

Big gcd(Big a, Big b) {
  if (a == 0) return b;
  if (b == 0) return a;
  // Count the powers of the common factor 2 in a and b; countr_zero gives the number of trailing zero bits
  int atimes = countr_zero(a);
  int btimes = countr_zero(b);
  int mintimes = min(atimes, btimes);
  a >>= atimes;
  for (;;) {
    // The common powers of 2 in a and b have been handled, so a cannot be even later
    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 code above follows the C++17 std::gcd implementations in libstdc++ and MSVC. For unsigned int and unsigned long long, Stein's algorithm is faster than the Euclidean algorithm when countr_zero can be computed extremely quickly, but it may be slower otherwise.

About countr_zero

GCC's built-in functions __builtin_ctz (32-bit) and __builtin_ctzll (64-bit) can replace countr_zero in the code above.
Starting with C++20, the <bit> header contains std::countr_zero.
If non-standard-library functions are unavailable and C++20 cannot be used, the following is an $O(1)$ implementation after preprocessing in the Word-RAM with multiplication model:

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 big-integer arithmetic, if the representation is similar to bitset, combining it with the implementation of countr_zero above gives time complexity O(n / w). If splitting by individual bits is inconvenient, the largest power of $2$ dividing the number must instead be found by brute force, and the complexity depends on the implementation. For example:

// A binary Big stored in little-endian order; every element must be enumerable
int countr_zero(Big a) {
  int ans = 0;
  for (auto x : a) {
    if (x != 0) {
      ans += 32;  // Bit width of each element type
    } else {
      return ans + countr_zero(x);
    }
  }
  return ans;
}

// Brute-force computation; if used, put it directly into gcd to reduce the constant factor
int countr_zero(Big a) {
  int ans = 0;
  while ((a & 1) == 0) {
    a >>= 1;
    ++ans;
  }
  return ans;
}

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

Greatest common divisor of multiple numbers

How do we find the greatest common divisor of multiple numbers? The answer is clearly a divisor of every number, and therefore also a divisor of every adjacent pair. By induction, removing two numbers, computing their gcd, and putting the result back does not change the desired 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 divisible by every integer in the set. 0 is a common multiple of every 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 by $\operatorname{lcm}(a,b)$ , or simply $[a,b]$ when there is no ambiguity.

For integers $a_1,\dots,a_n$ , their least common multiple is denoted by $\operatorname{lcm}(a_1,\dots,a_n)$ , or simply $[a_1,\dots,a_n]$ when there is no ambiguity.

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}}$ .

For $a$ and $b$ , we find that their greatest common divisor is

$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 is

$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 find the least common multiple of two numbers, first compute their greatest common divisor.

Multiple numbers

Once the $\gcd$ $O(1)$ $\gcd$ $\gcd$ of two numbers is known, their least common multiple can be computed in time. For multiple numbers, there is no need to compute a common gcd first: whenever the gcd of two numbers is computed in a sequence, put their least common multiple back into the sequence and continue with the remaining numbers.

Extended Euclidean algorithm

The extended Euclidean algorithm (EXGCD) is commonly used to find one 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 identity, $\gcd(a,b)=\gcd(b,a\bmod b)$ .

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

Also, $a\bmod b=a-(\lfloor\frac{a}{b}\rfloor\times b)$ .

Therefore $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)$

Comparing the coefficients of $a=a,b=b$ $x_1=y_2,y_1=x_2-\lfloor\frac{a}{b}\rfloor y_2$ and , we obtain and .

Substitute $x_2,y_2$ $b$ $0$ $x=1,y=0$ recursively until , where the recursion returns with .

Implementation

C++

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

Python

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 the $\gcd$ $x,y$ ; and are computed along the way.

Range analysis

The equation $ax+by=\gcd(a,b)$ $b\not= 0$ $|x|\le b,|y|\le a$ has infinitely many solutions, and some obviously exceed the range of the long long type. Fortunately, if , the solution produced by the extended Euclidean algorithm always satisfies and . We prove this property below.

Proof

When $\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$ , , so the recursion terminates at the next level. We obtain , and clearly .
When , assume . Since and , we have and ; this is at most . Therefore .

Writing the extended Euclidean algorithm iteratively

First, when $x = 1$ , $y = 0$ , $x_1 = 0$ , and $y_1 = 1$ , we clearly have:

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

which holds.

Given $a\bmod b = a - (\lfloor \frac{a}{b} \rfloor \times b)$ , let $q = \lfloor \frac{a}{b} \rfloor$ . Following the iterative gcd algorithm, one iteration can be written as:

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

Replace $a$ in the iteration by $ax + by = a$ and $b$ by $ax_1 + by_1 = b$ to obtain:

\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 gives the iterative extended Euclidean algorithm.

The iterative method avoids recursion, so it runs slightly faster than the recursive version.

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 inspect $a_1$ and $b_1$ , you will find that they have exactly the same values as in the iterative Euclidean algorithm. The identities $x \cdot a +y \cdot b =a_1$ and $x_1 \cdot a +y_1 \cdot b= b_1$ hold both before the while loop and after every iteration, so the algorithm necessarily computes the $\gcd$ correctly.

Finally, $a_1$ is the required $\gcd$ , and $x \cdot a +y \cdot b = g$ .

Matrix interpretation

For positive integers $a$ and $b$ , one Euclidean step, namely $\gcd(a,b)=\gcd(b,a\bmod b)$ , can be represented by the matrix

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

Here, the floor symbol $\lfloor c\rfloor$ denotes the greatest integer not exceeding $c$ . 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, giving

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

The condition $a\cdot x_1+b\cdot x_2=\gcd(a,b)$ is precisely the extended Euclidean algorithm. Multiplication by an identity matrix at the end does not change the result, suggesting that we can maintain a $2\times 2$ identity matrix from the start and write the following simpler iterative method:

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

Greatest common divisor

Definition

Euclidean algorithm

Process

Implementation

Properties

Subtractive Euclidean algorithm

Process

Stein algorithm optimization

Implementation

Greatest common divisor of multiple numbers

Least common multiple

Definition

Two numbers

Multiple numbers

Extended Euclidean algorithm

Process

Implementation

Range analysis

Writing the extended Euclidean algorithm iteratively

Matrix interpretation

Applications

References and links

Footnotes