Number Theory: GCD & Modular Arithmetic

Most of this course has measured algorithms against the size of their input: $n$ elements, $V$ vertices, $h$ levels of a tree. Number-theoretic algorithms break that habit. Their inputs are single integers, and the interesting cost is measured against the magnitude of those integers, or equivalently the number of bits needed to write them down. The oldest non-trivial algorithm we know, Euclid's, from around 300 BCE, belongs to this family, and it is still the right way to compute a greatest common divisor. This lesson develops it carefully, extends it to solve linear Diophantine equations, and uses both to lay the foundations of modular arithmetic, the arithmetic that underlies hashing, cryptography, and the math problems you will meet in practice.

Divisibility and the greatest common divisor

For integers $a$ and $d$ we say $d$ divides $a$ (written $d\mid a$) if $a=kd$ for some integer $k$. A common divisor of $a$ and $b$ is an integer dividing both. The greatest common divisor $\gcd (a,b)$ is the largest such integer, with the conventions $\gcd (a,0)=|a|$ and $\gcd (0,0)=0$. Throughout we take $a,b\ge 0$; signs only flip the answer's sign.¹

The naive way to compute $\gcd (a,b)$, factoring both numbers and multiplying the shared prime powers, is a mistake: integer factorization is believed to be hard, and the sieve and factorization methods that find those prime powers are themselves a separate study. Euclid's insight is that the factorization is never needed.

Euclid's algorithm

The entire algorithm rests on one recurrence.

The base case is immediate: every integer divides $0$, so the largest divisor of $a$ and $0$ is $a$ itself. The recurrence follows from a sharper claim: the two pairs share exactly the same set of common divisors, hence the same greatest one.

Because the second argument strictly shrinks ($a\mathrm{mod}b<b$) and stays non-negative, the recursion must terminate, and it terminates at a pair $(g,0)$ whose answer is $g=\gcd (a,b)$.

def gcd(first: int, second: int) -> int:
  """
    Greatest common divisor of `first` and `second`.\n
    Signs are ignored: the result is non-negative, with gcd(0, 0) = 0.\n
  """
  # signs don't affect divisors, so work with magnitudes.
  current: int = abs(first)
  divisor: int = abs(second)

  # replace (current, divisor) by (divisor, current mod divisor) until it's 0.
  while divisor != 0:
    current, divisor = divisor, current % divisor

  return current

def lcm(first: int, second: int) -> int:
  """
    Least common multiple, via lcm(a, b) = |a * b| / gcd(a, b).\n
    Defined as 0 when either argument is 0.\n
  """
  if first == 0 or second == 0:
    return 0

  # divide before multiplying to keep the intermediate value small.
  return abs(first // gcd(first, second) * second)

def gcd_of_all(values: list[int]) -> int:
  """
    Greatest common divisor of a whole list, folding gcd left to right.\n
    The gcd of an empty list is 0 (the identity for gcd).\n
  """
  result: int = 0
  for value in values:
    result = gcd(result, value)
  return result

Why it is fast

Each iteration replaces $(a,b)$ with $(b,a\mathrm{mod}b)$. The key fact is that two iterations at least halve the larger argument.

So after every two steps the first argument drops below half its value; the number of iterations is therefore $O(\log a)=O(\log \min (a,b))$ once the first swap orders the arguments.² Each iteration does one division on numbers of $O(\log a)$ bits, so the bit-complexity is polynomial in the input size, exponentially better than factoring. (For a refresher on this kind of logarithmic bound, see asymptotic analysis.)

The worst case is slow precisely because every quotient is the smallest it can be, $\lfloor a/b\rfloor =1$, so each step subtracts $b$ only once and the pair merely slides to the previous Fibonacci pair. A single larger quotient would collapse the chain far faster.

Each iteration simply replaces the pair by $(b,a\mathrm{mod}b)$ and recurses; tracing $\gcd (48,18)$ shows the second argument collapsing to $0$ in three steps.

Euclid, geometrically

Replacing $a\mathrm{mod}b$ by repeated subtraction gives the subtractive form of the algorithm, and it has a geometric reading: tile an $a\times b$ rectangle greedily with the largest squares that fit. Cut off a $b\times b$ square as many times as you can, then recurse on the leftover $b\times (a\mathrm{mod}b)$ strip. The side of the last square is the gcd.

Extended Euclid: Bézout's identity

Euclid tells us what the gcd is; the extended algorithm tells us how to build it out of $a$ and $b$.

So the back-substitution recurrence is

from typing import NamedTuple

class Bezout(NamedTuple):
  """
    A Bezout solution: gcd together with coefficients x, y such that\n
    a*x + b*y = gcd holds for the inputs (a, b) that produced it.\n
  """
  gcd: int
  x: int
  y: int

def extended_gcd(first: int, second: int) -> Bezout:
  """
    Greatest common divisor of `first` and `second` with Bezout coefficients\n
    x, y satisfying first*x + second*y == gcd. The returned gcd is\n
    non-negative; x and y may be negative.\n
    Implemented iteratively to avoid recursion depth on large inputs.\n
  """
  # track remainders alongside their bezout coefficient pairs.
  old_remainder, remainder = first, second
  old_x, current_x = 1, 0
  old_y, current_y = 0, 1

  # run euclid, carrying each coefficient through the same back-substitution.
  while remainder != 0:
    quotient: int = old_remainder // remainder
    old_remainder, remainder = remainder, old_remainder - quotient * remainder
    old_x, current_x = current_x, old_x - quotient * current_x
    old_y, current_y = current_y, old_y - quotient * current_y

  # inputs may leave the gcd negative; flip its sign and the coefficients.
  if old_remainder < 0:
    return Bezout(-old_remainder, -old_x, -old_y)

  return Bezout(old_remainder, old_x, old_y)

It performs the same divisions as plain Euclid, so it is also $O(\log \min (a,b))$. The table below traces $\gcd (240,46)$: the forward pass fills the remainder/quotient columns top-down, and the coefficients $(x,y)$ are filled bottom-up by the back-substitution recurrence, landing on Bézout coefficients for the original pair in the top row.

When does $ax+by=c$ have a solution?

Bézout characterizes exactly when the general linear Diophantine equation is solvable.

This is the predicate behind the Water and Jug Problem (can we measure $c$ liters using jugs of capacity $a$ and $b$? iff $g\mid c$ and $c\le a+b$) and Check if Point Is Reachable, where the reachable lattice is governed by the gcd of the allowed steps.

from typing import NamedTuple, Optional

from extended_gcd import extended_gcd

class DiophantineSolution(NamedTuple):
  """
    A particular solution (x, y) to a*x + b*y = c, plus the steps that\n
    generate every other solution: the general solution is\n
    (x + k*x_step, y + k*y_step) for every integer k.\n
  """
  x: int
  y: int
  x_step: int
  y_step: int

def solve_diophantine(
  coefficient_a: int, coefficient_b: int, target: int
) -> Optional[DiophantineSolution]:
  """
    Solve coefficient_a * x + coefficient_b * y = target in integers.\n
    Returns a `DiophantineSolution` describing the whole solution family, or\n
    None when no integer solution exists (i.e. gcd does not divide target).\n
    The degenerate all-zero-coefficient cases are handled explicitly.\n
  """
  # with both coefficients zero, only 0 = target (i.e. target 0) is solvable.
  if coefficient_a == 0 and coefficient_b == 0:
    if target != 0:
      return None
    return DiophantineSolution(0, 0, 1, 0)

  # bezout's gcd must divide the target for any integer solution to exist.
  bezout = extended_gcd(coefficient_a, coefficient_b)
  divisor: int = bezout.gcd
  if target % divisor != 0:
    return None

  # scale bezout's coefficients by target / gcd for a particular solution.
  scale: int = target // divisor
  particular_x: int = bezout.x * scale
  particular_y: int = bezout.y * scale

  # stepping along (b/g, -a/g) leaves a*x + b*y fixed, enumerating all others.
  x_step: int = coefficient_b // divisor
  y_step: int = -coefficient_a // divisor

  return DiophantineSolution(particular_x, particular_y, x_step, y_step)

def has_diophantine_solution(
  coefficient_a: int, coefficient_b: int, target: int
) -> bool:
  """
    Whether coefficient_a * x + coefficient_b * y = target is solvable in\n
    integers — exactly the predicate gcd(a, b) | target.\n
  """
  return solve_diophantine(coefficient_a, coefficient_b, target) is not None

Modular arithmetic

Fix a modulus $m>0$. We say $a$ is congruent to $b$ modulo $m$, written $a\equiv b(\mathrm{mod}m)\leftrightarrow m\mid (a-b),$ i.e. $a$ and $b$ leave the same remainder on division by $m$. Congruence is an equivalence relation, and it partitions the integers into $m$ residue classes$\{0,1,\dots ,m-1\}$. The decisive property is that the class operations are well-defined: if $a\equiv a^{\prime }$ and $b\equiv b^{\prime }$, then $a+b\equiv a^{\prime }+b^{\prime },a-b\equiv a^{\prime }-b^{\prime },a\cdot b\equiv a^{\prime }\cdot b^{\prime }(\mathrm{mod}m).$ So you may reduce mod $m$ at any point in a chain of $+$, $-$, $\times$ without changing the final residue, the foundation of every answer modulo $10^9+7$ problem and of combinatorics modulo a prime.³

When two coprime moduli are at play, the residue classes interlock perfectly: the pair $(x\mathrm{mod}3,x\mathrm{mod}5)$ pins $x$ down uniquely modulo $15$. The grid below tabulates that bijection, with $x=10(x\mathrm{mod}3)+6(x\mathrm{mod}5)$ landing in each cell, the constructive heart of the Chinese Remainder Theorem we revisit in combinatorics.

Modular inverse and linear congruences

A modular inverse of $a$ modulo $m$ is an integer $a^{-1}$ with $a\cdot a^{-1}\equiv 1(\mathrm{mod}m)$. It is what lets you divide by $a$.

There are two standard ways to produce the inverse:

1. Extended Euclid. Run $\text{Extended-Euclid}(a,m)$ to get $ax+my=1$. Reducing mod $m$ kills the $my$ term, leaving $ax\equiv 1(\mathrm{mod}m)$, so $x\mathrm{mod}m$ is the inverse. This works for any coprime modulus and costs $O(\log m)$.
2. Fermat's little theorem. When $m$ is prime, every $a\not\equiv 0$ is coprime to $m$, and $a^{m-1}\equiv 1(\mathrm{mod}m)$, hence $a^{-1}\equiv a^{m-2}(\mathrm{mod}m)$. Computed by fast exponentiation in $O(\log m)$ multiplications, the subject of the next lesson, on modular exponentiation and primality.

The reason an inverse exists exactly when $\gcd (a,m)=1$ is visible directly: multiplying every nonzero residue by such an $a$ permutes them, so some residue must land on $1$, and that residue is $a^{-1}$. Below, multiplying $\{1,\dots ,6\}$ by $3$ modulo $7$ shuffles the set, and the arrow into $1$ comes from $5$, so $3^{-1}\equiv 5$.

from typing import Optional

from extended_gcd import extended_gcd

def mod_inverse(value: int, modulus: int) -> Optional[int]:
  """
    The inverse of `value` modulo `modulus`, as a representative in\n
    [0, modulus), or None when `value` is not coprime to `modulus`.\n
    `modulus` must be at least 1; the only unit modulo 1 is 0.\n
  """
  if modulus <= 0:
    raise ValueError("modulus must be positive")

  # everything is congruent to 0 modulo 1, and 0 is its own inverse there.
  if modulus == 1:
    return 0

  # an inverse exists only when value is coprime to the modulus.
  bezout = extended_gcd(value % modulus, modulus)
  if bezout.gcd != 1:
    return None

  # normalize x into [0, modulus); bezout.x can be negative.
  return ((bezout.x % modulus) + modulus) % modulus

The same machinery solves the general linear congruence $ax\equiv b(\mathrm{mod}m)$. Let $g=\gcd (a,m)$.

When $g=1$ this reduces to multiply both sides by $a^{-1}$ and yields the single solution $x\equiv a^{-1}b(\mathrm{mod}m)$, the everyday case. The general count $g$ applies when the modulus and coefficient share a factor.

from extended_gcd import extended_gcd
from mod_inverse import mod_inverse

def solve_linear_congruence(
  coefficient: int, target: int, modulus: int
) -> list[int]:
  """
    All solutions to coefficient * x ≡ target (mod modulus), as a sorted list\n
    of representatives in [0, modulus). Returns an empty list when there is no\n
    solution (i.e. gcd(coefficient, modulus) does not divide target).\n
    `modulus` must be at least 1.\n
  """
  if modulus <= 0:
    raise ValueError("modulus must be positive")

  # solvable only when gcd(coefficient, modulus) divides the target.
  divisor: int = extended_gcd(coefficient, modulus).gcd
  if target % divisor != 0:
    return []

  # divide through by the gcd to a coprime congruence mod m/g.
  reduced_modulus: int = modulus // divisor
  reduced_coefficient: int = (coefficient // divisor) % reduced_modulus
  reduced_target: int = (target // divisor) % reduced_modulus

  # the reduced coefficient is coprime to m/g, so its inverse always exists.
  inverse = mod_inverse(reduced_coefficient, reduced_modulus)
  assert inverse is not None  # coprimality guarantees the inverse exists.
  base_solution: int = (reduced_target * inverse) % reduced_modulus

  # the g distinct solutions are spaced m/g apart; return them sorted.
  solutions: list[int] = [
    (base_solution + step * reduced_modulus) % modulus
    for step in range(divisor)
  ]
  solutions.sort()

  return solutions

Worked example (a linear congruence with a shared factor). Solve $6x\equiv 8(\mathrm{mod}14)$. Here $a=6$, $b=8$, $m=14$, and $g=\gcd (6,14)=2$. Since $g=2$ divides $b=8$, the congruence is solvable, and the theorem promises exactly $g=2$ solutions modulo $14$. Divide the whole congruence through by $g$: with $a^{\prime }=3$, $b^{\prime }=4$, $m^{\prime }=7$ we solve $3x\equiv 4(\mathrm{mod}7)$. The inverse $3^{-1}\equiv 5(\mathrm{mod}7)$ (from the permutation figure above), so $x_0\equiv 5\cdot 4\equiv 20\equiv 6(\mathrm{mod}7)$. Lifting back to modulus $14$, the two solutions are $x_0$ and $x_0+m^{\prime }=6$ and $13$. Checking: $6\cdot 6=36\equiv 8(\mathrm{mod}14)$ and $6\cdot 13=78\equiv 8(\mathrm{mod}14)$, both correct.

GCD refinements and where Bézout matters

Euclid's algorithm is the oldest non-trivial algorithm in continuous use, and the modern refinements of it are worth knowing.

Binary GCD (Stein's algorithm). Division is expensive on hardware that lacks a fast divide, and each Euclid step needs a modulo. In 1967 Josef Stein published a variant that uses only subtraction, comparison, and shifts — no division at all.⁴ It rests on three facts: $\gcd (2a,2b)=2\gcd (a,b)$ (pull out a common factor of two), $\gcd (2a,b)=\gcd (a,b)$ when $b$ is odd (a factor of two in one argument alone is irrelevant to an odd gcd), and $\gcd (a,b)=\gcd (|a-b|,\min (a,b))$ for $a,b$ odd. Stripping factors of two is a single shift instruction, so binary GCD is often faster in practice than Euclid despite touching the same $O(\log )$ number of bits. On $\gcd (24,18)$: pull out one common $2$ to reach $2\gcd (12,9)$; $12$ is even and $9$ odd, so drop that factor to $2\gcd (6,9)$, then $2\gcd (3,9)$; both odd now, subtract to $2\gcd (6,3)=2\gcd (3,3)=2\cdot 3=6$.

Bit complexity, done honestly. Counting each division as one step gives $O(\log \min (a,b))$ operations, but on numbers of $n=O(\log a)$ bits, one schoolbook division already costs $O(n^2)$ bit operations, so Euclid is $O(n^3)$ bit operations in the naive accounting. The half-GCD algorithm, using fast multiplication, computes a gcd in $O(n{\log }^{2}n\log \log n)$ bit operations by a Knuth–Schönhage divide-and-conquer that processes the high-order bits in one batch,⁵ the same asymptotic class as multiplication itself. This is what large-integer libraries (GMP) actually run for big inputs.

Where Bézout coefficients matter. The extended algorithm's coefficients $(x,y)$ serve directly as the modular inverse ($x\mathrm{mod}m$), the CRT reconstruction weights, and the private exponent in RSA key generation ($d\equiv e^{-1}(\mathrm{mod}\phi (N))$ is one extended-Euclid call). Every time a cryptographic library inverts modulo the group order, it is running the algorithm on this page.

Takeaways

• $\gcd (a,b)$ is computed by Euclid's algorithm via $\gcd (a,b)=\gcd (b,a\mathrm{mod}b)$, base $\gcd (a,0)=a$; the recurrence is exact because $(a,b)$ and $(b,a\mathrm{mod}b)$ have identical common divisors.
• Euclid runs in $O(\log \min (a,b))$ iterations because two steps at least halve the argument; the worst case is consecutive Fibonacci numbers.
• Extended Euclid returns Bézout coefficients $x,y$ with $ax+by=\gcd (a,b)$, and $ax+by=c$ is solvable iff $\gcd (a,b)\mid c$, the test behind Water-and-Jug and Check-if-Point-Is-Reachable.
• Modular arithmetic is arithmetic on residue classes; $+,-,\times$ are well-defined, but division requires an inverse and overflow must be guarded.
• A modular inverse $a^{-1}\mathrm{mod}m$ exists iff $\gcd (a,m)=1$, found by extended Euclid, or by Fermat ($a^{m-2}$) when $m$ is prime; the linear congruence $ax\equiv b(\mathrm{mod}m)$ is solvable iff $\gcd (a,m)\mid b$, with exactly $\gcd (a,m)$ solutions.