Euclid’s algorithm in Python: use remainders instead of repeated subtraction

Problem. A greatest-common-divisor helper can waste work by subtracting the smaller number one copy at a time, even when one remainder would discard the whole stack. Baseline complexity. The literal subtraction loop is linear in the value of the larger operand in its bad cases: for (600001, 1), it performs 600000 loop iterations. Solution. Euclid’s algorithm replaces a pair with its divisor and remainder, preserving the shared divisors while making much larger jumps. Measured result. On this CPython 3.13.5 aarch64 run, (600001, 1) takes 600000 subtraction updates but one remainder update. In a matched 1000-call benchmark of (6001, 1), the subtraction version’s best per-call result was 0.001259863 seconds and the remainder version’s was 0.000000965 seconds, a 1305.5× ratio.

Reproduce the result

Complete code

from __future__ import annotations
import json
import math
import platform
import sys
import timeit
def require_nonnegative_integers(a: int, b: int) -> None:
    if isinstance(a, bool) or isinstance(b, bool) or not isinstance(a, int) or not isinstance(b, int):
        raise TypeError("a and b must be integers, excluding bool")
    if a < 0 or b < 0:
        raise ValueError("a and b must be non-negative")
    if a == 0 and b == 0:
        raise ValueError("gcd(0, 0) is not defined by this example")
def gcd_by_subtraction(a: int, b: int, *, count_steps: bool = False) -> int | tuple[int, int]:
    """A deliberately literal baseline: remove one smaller operand at a time."""
    require_nonnegative_integers(a, b)
    steps = 0
    if a == 0:
        return (b, steps) if count_steps else b
    if b == 0:
        return (a, steps) if count_steps else a
    while a != b:
        if a > b:
            a -= b
        else:
            b -= a
        steps += 1
    return (a, steps) if count_steps else a
def gcd_by_remainder(a: int, b: int, *, count_steps: bool = False) -> int | tuple[int, int]:
    """Euclid's algorithm: replace a pair with its divisor and remainder."""
    require_nonnegative_integers(a, b)
    steps = 0
    while b:
        a, b = b, a % b
        steps += 1
    return (a, steps) if count_steps else a
def expect_error(function, arguments: tuple[object, object], error_type: type[Exception]) -> bool:
    try:
        function(*arguments)
    except error_type:
        return True
    return False
def main() -> None:
    correctness_cases = [(48, 18), (21, 21), (17, 5), (0, 9), (9, 0), (600_001, 1)]
    correctness = []
    for a, b in correctness_cases:
        subtraction_value, subtraction_steps = gcd_by_subtraction(a, b, count_steps=True)
        remainder_value, remainder_steps = gcd_by_remainder(a, b, count_steps=True)
        reference = math.gcd(a, b)
        assert subtraction_value == remainder_value == reference
        correctness.append(
            {
                "pair": [a, b],
                "gcd": reference,
                "subtraction_steps": subtraction_steps,
                "remainder_steps": remainder_steps,
            }
        )
    edge_checks = {
        "negative_rejected": expect_error(gcd_by_remainder, (-1, 7), ValueError),
        "bool_rejected": expect_error(gcd_by_remainder, (True, 7), TypeError),
        "both_zero_rejected": expect_error(gcd_by_remainder, (0, 0), ValueError),
    }
    assert all(edge_checks.values())
    benchmark_pair = (6_001, 1)
    benchmark_batch_calls = 1_000
    baseline_batch_samples = timeit.repeat(
        lambda: gcd_by_subtraction(*benchmark_pair),
        number=benchmark_batch_calls,
        repeat=5,
    )
    remainder_batch_samples = timeit.repeat(
        lambda: gcd_by_remainder(*benchmark_pair),
        number=benchmark_batch_calls,
        repeat=5,
    )
    baseline_per_call_samples = [sample / benchmark_batch_calls for sample in baseline_batch_samples]
    remainder_per_call_samples = [sample / benchmark_batch_calls for sample in remainder_batch_samples]
    ratio = min(baseline_per_call_samples) / min(remainder_per_call_samples)
    report = {
        "benchmark_pair": list(benchmark_pair),
        "correctness_cases": correctness,
        "edge_checks": edge_checks,
        "benchmark_batch_calls": benchmark_batch_calls,
        "baseline_batch_samples_seconds": baseline_batch_samples,
        "baseline_per_call_samples_seconds": baseline_per_call_samples,
        "remainder_batch_samples_seconds": remainder_batch_samples,
        "remainder_per_call_samples_seconds": remainder_per_call_samples,
        "best_baseline_over_remainder_ratio": ratio,
        "python": sys.version,
        "platform": platform.platform(),
        "implementation": platform.python_implementation(),
        "timer": "timeit.repeat callable; five samples; equal 1000-call batches for both functions; GC disabled by timeit during timing",
    }
    print(json.dumps(report, indent=2, sort_keys=True))
if __name__ == "__main__":
    main()

Output

{
  "baseline_batch_samples_seconds": [
    1.2598634539172053,
    1.2717048190534115,
    1.298604816896841,
    1.3680981991346925,
    1.3163886698894203
  ],
  "baseline_per_call_samples_seconds": [
    0.0012598634539172054,
    0.0012717048190534115,
    0.0012986048168968408,
    0.0013680981991346926,
    0.0013163886698894203
  ],
  "benchmark_batch_calls": 1000,
  "benchmark_pair": [
    6001,
    1
  ],
  "best_baseline_over_remainder_ratio": 1305.4686685676252,
  "correctness_cases": [
    {
      "gcd": 6,
      "pair": [
        48,
        18
      ],
      "remainder_steps": 3,
      "subtraction_steps": 4
    },
    {
      "gcd": 21,
      "pair": [
        21,
        21
      ],
      "remainder_steps": 1,
      "subtraction_steps": 0
    },
    {
      "gcd": 1,
      "pair": [
        17,
        5
      ],
      "remainder_steps": 3,
      "subtraction_steps": 6
    },
    {
      "gcd": 9,
      "pair": [
        0,
        9
      ],
      "remainder_steps": 1,
      "subtraction_steps": 0
    },
    {
      "gcd": 9,
      "pair": [
        9,
        0
      ],
      "remainder_steps": 0,
      "subtraction_steps": 0
    },
    {
      "gcd": 1,
      "pair": [
        600001,
        1
      ],
      "remainder_steps": 1,
      "subtraction_steps": 600000
    }
  ],
  "edge_checks": {
    "bool_rejected": true,
    "both_zero_rejected": true,
    "negative_rejected": true
  },
  "implementation": "CPython",
  "platform": "Linux-6.12.47+rpt-rpi-v8-aarch64-with-glibc2.41",
  "python": "3.13.5 (main, Jul 15 2026, 20:25:40) [GCC 14.2.0]",
  "remainder_batch_samples_seconds": [
    0.0010199351236224174,
    0.0009650660213083029,
    0.0010473981965333223,
    0.001026269979774952,
    0.0010280651040375233
  ],
  "remainder_per_call_samples_seconds": [
    1.0199351236224174e-06,
    9.650660213083028e-07,
    1.0473981965333224e-06,
    1.026269979774952e-06,
    1.0280651040375233e-06
  ],
  "timer": "timeit.repeat callable; five samples; equal 1000-call batches for both functions; GC disabled by timeit during timing"
}

Environment

The executable recorded CPython 3.13.5 on Linux 6.12.47+rpt-rpi-v8, aarch64, with no third-party dependencies.

The reference check uses math.gcd, which Python documents as the greatest common divisor of its integer arguments.^[1]

Methodology

The program checks both implementations against math.gcd on six non-negative input pairs before timing, and it checks rejection of negatives, booleans and (0, 0).

It then times (6001, 1) five times with timeit.repeat, using an equal 1000-call batch for each function and dividing each batch back to a per-call value.

Python describes timeit as a way to time small pieces of code, and documents that it temporarily disables garbage collection during timing by default.^[2]

The complete raw vectors, operation counts, timer policy and platform string are in Output.

The problem in practical code

GCD and LCM methods also appear in application problems involving equal grouping and common multiples.^[5]

Python’s standard library already supplies math.gcd; the point of the two functions here is to make the work visible, rather than to replace that library function.^[1]

The subtraction version is a useful first sketch because it keeps the invariant easy to see.

If a > b, replacing a with a - b leaves the set of common divisors unchanged.

Its awkward habit appears when b is small.

With (600001, 1), each pass removes one. The loop does that 600000 times before the operands meet.

Baseline complexity

For positive integers, this literal baseline performs one subtraction per loop and can take a number of loops proportional to the larger numeric value.

The pair (n, 1) demonstrates the boundary exactly: it needs n - 1 subtraction updates.

That is a claim about this implementation’s loop count, not a claim that every GCD routine costs that much memory or CPU time.

The function holds only its two working integers and a step counter, so its algorithmic auxiliary state is constant for ordinary fixed-width values.

Python integers have variable-size representations, and the article does not use this as a process-memory measurement.

Remainders make the same invariant do more work

The Euclidean algorithm finds the GCD by repeatedly taking remainders of consecutive terms.^[3]

For a = bq + r, a number that divides both a and b also divides r, so gcd(a, b) equals gcd(b, r).^[3]

One modulo operation can therefore replace many individual subtractions.

Lamé’s theorem says that the number of Euclidean division steps never exceeds five times the number of decimal digits in the smaller positive operand.^[4]

That familiar logarithmic step bound counts remainder updates. It does not turn arbitrary-precision division into a free operation, and it says nothing about a particular Python build’s machine instructions.

Comparison diagram for GCD of 600001 and 1: repeated subtraction takes 600000 loop iterations; Euclid's remainder method takes one; both return 1. — Same answer, fewer laps. The diagram records loop counts; the executable keeps the stopwatch honest.

What the measurement shows

Both functions returned the same GCD as math.gcd for all six checked pairs, including zeros in either position.

For (6001, 1), the equal-batch trace reports 6000 subtraction updates per call and one remainder update.

The best subtraction sample was 0.001259863 seconds per call. The best remainder sample was 0.965 microseconds per call after normalising the same 1000-call batch, giving the recorded 1305.5× ratio.

That ratio belongs to this host, this Python build, this input pair and this equal-batch timer arrangement.

The pair was selected to expose the subtraction loop’s boundary, so it should not be read as an average speed-up for routine inputs.

The raw five-sample vectors matter here. A single neat-looking number would have been a slightly suspicious amount of tidiness.

When the technique is not the decision

Use math.gcd in production Python when a standard GCD is what the program needs.^[1]

The two teaching functions reject negative inputs and (0, 0) by design, while math.gcd has its own documented interface and accepts multiple integer arguments.^[1]

Repeated subtraction can still be adequate for a hand-worked example, very small values, or an environment where division is unavailable.

Euclid’s method also does not solve every number-theory task by itself. Extended Euclid is used when coefficients are required, because it returns integers that express the GCD as a linear combination of the inputs.^[3]^[6] A batch workload may need a different representation or a library routine.