The Euclidean Algorithm and GCD Computation

Executive summary

Three algorithms for one operation, chosen by operand size

The greatest common divisor is computed by repeated remaindering, and the number of steps is bounded by roughly 2.078 log₁₀ of the smaller operand — a bound attained exactly by consecutive Fibonacci numbers. For single-word operands the binary algorithm replaces division with shifts and subtractions; for genuinely multiprecision operands Lehmer's method performs many steps using only single-precision arithmetic before touching the full operands.

Learning objectives

State Euclid's algorithm and prove why it terminates.
Identify the worst-case input family and the resulting step bound.
Explain how the binary GCD avoids division entirely.
Describe how Lehmer's method batches steps into single-precision work.
Select the appropriate GCD variant for a given operand size.

Section 01The classical algorithm

The algorithm rests on one identity: for b ≠ 0,

gcd(a, b) = gcd(b, a mod b)

The second argument strictly decreases and stays non-negative, so the process terminates. The last non-zero remainder is the GCD.

AlgorithmEuclid's algorithmin: a, b ∈ ℤ → out: gcd(a, b) ≥ 0

If b = 0, return |a|.
Set r ← a mod b. The single expensive operation per iteration.
Set a ← b and b ← r.
Return to step 1.

Why the remainder, not repeated subtraction

Replacing a mod b with repeated subtraction is correct but can take Θ(a/b) steps in a single iteration. The division is what makes the algorithm logarithmic rather than linear.

Section 02Worst case and average behaviour

Lamé's theorem gives the worst case: the number of division steps is at most about five times the number of decimal digits of the smaller operand. The bound is tight, and it is attained by consecutive Fibonacci numbers — the inputs for which every quotient is 1 and therefore every step removes the least possible information.

O(log min(a,b))division steps

≈ 2.078 log₁₀Lamé bound per decimal digit

O(n²)bit operations, classical

Average behaviour is better than the worst case and is governed by the Gauss–Kuzmin distribution of continued fraction quotients: most quotients are small, but the mean number of steps is still logarithmic in the input. The practical consequence is that GCD cost is dominated by the cost of the divisions, not by their count.

GCD and continued fractions are the same computation

The sequence of quotients produced by Euclid's algorithm on (a, b) is exactly the continued fraction expansion of a/b. Every statement about one transfers to the other, which is why the worst case involves Fibonacci numbers — the ratio of consecutive Fibonacci numbers has the slowest-converging continued fraction there is.

Section 03The binary GCD

Stein's algorithm eliminates division. It uses three facts: the GCD of two even numbers is twice the GCD of their halves; a factor of 2 can be removed from an even operand when the other is odd; and the GCD of two odd numbers equals the GCD of the smaller and half their difference.

AlgorithmBinary GCD (Stein)in: a, b ≥ 0 → out: gcd(a, b)

If a = 0 return b; if b = 0 return a.
Let k be the number of common factors of 2; set a ← a/2^k, b ← b/2^k. One count-trailing-zeros instruction each.
While a is even, set a ← a/2.
Loop: while b is even, set b ← b/2.
If a > b, swap a and b.
Set b ← b − a. Both odd, so the difference is even.
If b = 0, return a · 2^k; otherwise repeat from step 4.

The binary algorithm wins on single-word and short multiword operands, and on any architecture where division is markedly slower than shifting. It loses on very large operands, where the number of iterations — proportional to the total bit length — dominates.

Section 04Lehmer's method for large operands

On multiprecision operands both classical and binary GCD spend nearly all their time touching full-length numbers to extract a small amount of information. Lehmer's observation is that the sequence of quotients is usually determined by the leading limbs alone.

Stage 01Extract leading wordsTake the top limb (or two) of each operand as single-precision approximations.
Stage 02Simulate stepsRun Euclid on the approximations, tracking a 2×2 transformation matrix, while the quotients from the upper and lower bounds agree.
Stage 03Apply in bulkMultiply the full-length operands by the accumulated matrix — one bulk update for many simulated steps.
Stage 04Fall back when ambiguousIf the bounds disagree at the first step, perform one exact multiprecision division and resume.

Why this is a real speedup

Each bulk update performs a linear number of limb operations but advances the computation by as many steps as were simulated in single precision. The asymptotic complexity class is unchanged; the constant factor improves substantially, which is what matters in practice.

Section 05Selection guidance

How large are the operands?
- Single word Binary GCD — no division unit needed, minimal overhead.
- A few words Binary GCD or classical Euclid — measure; the crossover is hardware-dependent.
- Many words Lehmer — batching pays for itself well before a hundred limbs.
- Very large, and asymptotics matter Half-GCD / subquadratic methods — divide-and-conquer variants that approach the cost of multiplication.
Are Bézout coefficients needed?
- Yes Extended Euclidean algorithm — see the dedicated page; the extension changes the cost model.
- No Plain GCD — do not compute cofactors you will discard.

Do not compute more than you need

Requesting Bézout coefficients when only the GCD is wanted roughly doubles the work and, for multiprecision inputs, adds substantial memory traffic. This is one of the most common avoidable inefficiencies in number-theoretic code.

ReferenceFrequently asked questions

Is gcd(0, 0) defined?

By the standard convention it is 0, which keeps the identity gcd(a, 0) = |a| valid without a special case. Implementations should document their choice, because the alternative convention of raising an error also appears.

Why is the Fibonacci sequence the worst case?

Because every quotient in the expansion is 1, so each step reduces the operands by the smallest possible amount. Any input requiring the maximum number of steps for its size must have all quotients equal to 1, which forces the Fibonacci pattern.

Does the binary GCD extend to Bezout coefficients?

Yes — the extended binary GCD exists — but the bookkeeping is more intricate than in the classical extended algorithm and the halving steps require care with the cofactors. Most libraries use the classical extended algorithm, or Lehmer's method extended, for that purpose.

NavigateContinue in this stream

Curated next steps from this page. The site also surfaces algorithmically related reading below.

ProvenanceSources and further reading

Knuth volume 2 gives the definitive analysis, including the average-case behaviour and Lehmer's method. Cohen's GTM 138 treats the multiprecision variants in an algorithmic style.

This page is an original KEVOS explanatory article. It presents the underlying mathematics — definitions, algorithms, complexity results and selection criteria — in KEVOS editorial voice. No text is reproduced from any copyrighted source. Where numerical tables are relevant, KEVOS links to live authoritative databases rather than republishing static values.