826: GCD, LCM, and the Euclidean Algorithm

Difficulty: 3 Level: Intermediate The oldest non-trivial algorithm (300 BC): compute GCD in O(log min(a,b)) — the cornerstone of all modular arithmetic and number theory.

The Problem This Solves

GCD and LCM appear in every corner of computer science: fraction simplification, scheduling (when do two periodic events coincide?), modular inverse computation, Chinese Remainder Theorem, RSA key generation, and competitive programming problem after problem. The Euclidean algorithm is the prerequisite for all of these. The seemingly simple `gcd(a, b) = gcd(b, a mod b)` recurrence is deceptively powerful: it runs in O(log min(a, b)) steps because the remainder shrinks by at least half every two steps (Fibonacci analysis). For 64-bit integers, this means at most ~93 iterations. No algorithm with fewer comparisons exists. The Extended Euclidean algorithm (see #832) produces Bézout coefficients — integers x, y such that ax + by = gcd(a, b) — enabling modular inverses and CRT. Understanding plain GCD first makes the extension straightforward.

The Intuition

`gcd(a, b)` = any common divisor of a and b also divides `a mod b` (because `a mod b = a - ⌊a/b⌋ × b`), and the common divisors of (b, a mod b) are exactly the common divisors of (a, b). So `gcd(a, b) = gcd(b, a mod b)`. Base case: `gcd(a, 0) = a`. Efficiency: `a mod b < a/2` whenever `b ≤ a/2`, so the problem halves every two recursive calls. LCM: `lcm(a, b) = a / gcd(a, b) × b`. Divide before multiplying to prevent overflow — a common bug to know and avoid.

How It Works in Rust

// Recursive: mirrors OCaml's one-liner elegantly
fn gcd(a: u64, b: u64) -> u64 {
 if b == 0 { a } else { gcd(b, a % b) }
}

// Iterative: preferred for large inputs (no stack concerns)
// Rust does NOT guarantee tail-call optimization — always provide iterative version
fn gcd_iter(mut a: u64, mut b: u64) -> u64 {
 while b != 0 {
     let r = a % b;  // Save remainder before overwriting
     a = b;
     b = r;
 }
 a
}

// LCM: divide first to prevent overflow — critical for large u64 values
fn lcm(a: u64, b: u64) -> u64 {
 if a == 0 || b == 0 { return 0; }
 a / gcd(a, b) * b  // NOT: a * b / gcd(a, b) — would overflow for large a, b
}

// GCD of a slice: fold with neutral element 0 (gcd(0, x) = x)
fn gcd_slice(xs: &[u64]) -> u64 {
 xs.iter().fold(0u64, |acc, &x| gcd(acc, x))
}

// Binary GCD (Stein's): replaces modulo with bit shifts
// Faster on CPUs where division is expensive; same O(log n) asymptotically
fn binary_gcd(mut a: u64, mut b: u64) -> u64 {
 if a == 0 { return b; }
 if b == 0 { return a; }
 let shift = (a | b).trailing_zeros();
 a >>= a.trailing_zeros();
 loop {
     b >>= b.trailing_zeros();          // Remove factors of 2
     if a > b { std::mem::swap(&mut a, &mut b); }
     b -= a;                             // b = b - a, now even
     if b == 0 { break; }
 }
 a << shift
}

Note: Rust does not guarantee tail-call optimization. The recursive `gcd` is fine in practice (max ~93 stack frames for u64), but for production code handling untrusted inputs, prefer `gcd_iter`.

What This Unlocks

Modular inverse and cryptography: GCD = 1 is the prerequisite for modular inverse; the Extended Euclidean algorithm (built on this) is the foundation of RSA, ElGamal, and ECDSA.
Fraction arithmetic: Simplify `a/b` to lowest terms with `gcd(a, b)` in O(log min(a,b)) — used in symbolic math, exact arithmetic, and rendering pipelines.
Scheduling and LCM: "When do two tasks with period p1 and p2 synchronize?" → LCM(p1, p2). Used in real-time OS scheduling and animation frame synchronization.

Key Differences

Concept	OCaml	Rust
Recursive GCD	`let rec gcd a b = if b=0 then a else gcd b (a mod b)`	Identical structure; same risk of deep recursion
Tail recursion	OCaml compiler optimizes tail calls	Rust does NOT guarantee TCO; use iterative for safety
LCM overflow	`a b / gcd a b` (risky)	`a / gcd(a, b) b` — divide first
Slice fold	`List.fold_left gcd 0`	`iter().fold(0u64, gcd)`
Binary GCD	External implementation	`trailing_zeros()` makes Stein's idiomatic

/// GCD, LCM, and the Euclidean Algorithm.
///
/// gcd(a, b) = gcd(b, a % b) until b = 0.
/// Converges in O(log min(a, b)) steps.

/// Recursive GCD — mirrors OCaml's one-liner elegantly.
fn gcd(a: u64, b: u64) -> u64 {
    if b == 0 { a } else { gcd(b, a % b) }
}

/// Iterative GCD — no stack concerns for large inputs.
fn gcd_iter(mut a: u64, mut b: u64) -> u64 {
    while b != 0 {
        let r = a % b;
        a = b;
        b = r;
    }
    a
}

/// LCM: divide first to prevent overflow.
fn lcm(a: u64, b: u64) -> u64 {
    if a == 0 || b == 0 { return 0; }
    a / gcd(a, b) * b
}

/// GCD of a slice — fold with neutral element 0.
fn gcd_slice(xs: &[u64]) -> u64 {
    xs.iter().fold(0u64, |acc, &x| gcd(acc, x))
}

/// LCM of a slice — fold with neutral element 1.
fn lcm_slice(xs: &[u64]) -> u64 {
    xs.iter().fold(1u64, |acc, &x| lcm(acc, x))
}

/// Binary GCD (Stein's algorithm) — replaces modulo with bit shifts.
fn binary_gcd(mut a: u64, mut b: u64) -> u64 {
    if a == 0 { return b; }
    if b == 0 { return a; }
    let shift = (a | b).trailing_zeros();
    a >>= a.trailing_zeros();
    loop {
        b >>= b.trailing_zeros();
        if a > b { std::mem::swap(&mut a, &mut b); }
        b -= a;
        if b == 0 { break; }
    }
    a << shift
}

/// Trace the Euclidean algorithm steps.
fn gcd_trace(mut a: u64, mut b: u64) -> u64 {
    println!("Tracing gcd({a}, {b}):");
    while b != 0 {
        let r = a % b;
        println!("  gcd({a}, {b}) → remainder {r}");
        a = b;
        b = r;
    }
    println!("  = {a}");
    a
}

fn main() {
    gcd_trace(48, 18);

    println!("\ngcd(48, 18) = {} (expected 6)", gcd(48, 18));
    println!("lcm(4, 6)   = {} (expected 12)", lcm(4, 6));
    println!("lcm(12, 18) = {} (expected 36)", lcm(12, 18));
    println!("gcd_slice([12,18,24]) = {} (expected 6)", gcd_slice(&[12, 18, 24]));
    println!("lcm_slice([2,3,4])    = {} (expected 12)", lcm_slice(&[2, 3, 4]));
    println!("binary_gcd(48, 18) = {} (expected 6)", binary_gcd(48, 18));
    println!("\nLCM of 1..10 = {}", lcm_slice(&[1,2,3,4,5,6,7,8,9,10])); // 2520
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_gcd_basic() {
        assert_eq!(gcd(48, 18), 6);
        assert_eq!(gcd(18, 48), 6); // Commutative
    }

    #[test]
    fn test_gcd_coprime() {
        assert_eq!(gcd(13, 7), 1);
    }

    #[test]
    fn test_gcd_same() {
        assert_eq!(gcd(12, 12), 12);
    }

    #[test]
    fn test_gcd_zero() {
        assert_eq!(gcd(0, 5), 5);
        assert_eq!(gcd(5, 0), 5);
    }

    #[test]
    fn test_lcm_basic() {
        assert_eq!(lcm(4, 6), 12);
        assert_eq!(lcm(12, 18), 36);
    }

    #[test]
    fn test_lcm_coprime() {
        assert_eq!(lcm(3, 7), 21);
    }

    #[test]
    fn test_gcd_iter_matches() {
        for a in 0..50u64 {
            for b in 0..50u64 {
                assert_eq!(gcd(a, b), gcd_iter(a, b));
            }
        }
    }

    #[test]
    fn test_binary_gcd_matches() {
        for a in 0..50u64 {
            for b in 0..50u64 {
                assert_eq!(gcd(a, b), binary_gcd(a, b),
                    "binary_gcd({a}, {b}) mismatch");
            }
        }
    }

    #[test]
    fn test_gcd_slice() {
        assert_eq!(gcd_slice(&[12, 18, 24]), 6);
        assert_eq!(gcd_slice(&[7, 14, 21]), 7);
    }

    #[test]
    fn test_lcm_slice() {
        assert_eq!(lcm_slice(&[2, 3, 4]), 12);
        assert_eq!(lcm_slice(&[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]), 2520);
    }
}

(* GCD, LCM, and Euclidean Algorithm in OCaml *)

(* The classic one-liner — functional elegance at its finest *)
let rec gcd (a : int) (b : int) : int =
  if b = 0 then a else gcd b (a mod b)

(* LCM via GCD — divide first to avoid overflow *)
let lcm (a : int) (b : int) : int =
  if a = 0 || b = 0 then 0
  else (a / gcd a b) * b

(* GCD of a list — neutral element is 0 since gcd(a,0)=a *)
let list_gcd (xs : int list) : int =
  List.fold_left gcd 0 xs

let list_lcm (xs : int list) : int =
  List.fold_left lcm 1 xs

(* Trace the steps of Euclidean algorithm *)
let gcd_trace (a : int) (b : int) : unit =
  Printf.printf "gcd(%d, %d):\n" a b;
  let a = ref a and b = ref b in
  while !b <> 0 do
    Printf.printf "  gcd(%d, %d) → rem = %d\n" !a !b (!a mod !b);
    let r = !a mod !b in
    a := !b;
    b := r
  done;
  Printf.printf "  = %d\n" !a

(* Binary GCD (Stein's algorithm) — avoids modulo *)
let rec binary_gcd (a : int) (b : int) : int =
  match (a, b) with
  | (0, b) -> b
  | (a, 0) -> a
  | (a, b) when a mod 2 = 0 && b mod 2 = 0 -> 2 * binary_gcd (a / 2) (b / 2)
  | (a, b) when a mod 2 = 0 -> binary_gcd (a / 2) b
  | (a, b) when b mod 2 = 0 -> binary_gcd a (b / 2)
  | (a, b) -> binary_gcd (abs (a - b)) (min a b)

let () =
  gcd_trace 48 18;
  Printf.printf "\ngcd(48, 18) = %d  (expected 6)\n" (gcd 48 18);
  Printf.printf "lcm(4, 6)   = %d  (expected 12)\n" (lcm 4 6);
  Printf.printf "lcm(12, 18) = %d  (expected 36)\n" (lcm 12 18);
  Printf.printf "gcd [12;18;24] = %d  (expected 6)\n" (list_gcd [12;18;24]);
  Printf.printf "lcm [2;3;4]    = %d  (expected 12)\n" (list_lcm [2;3;4]);
  Printf.printf "binary_gcd(48, 18) = %d\n" (binary_gcd 48 18)