Lab 07 — Extended Euclid, Modular Inverse & CRT

Goal

Implement extended gcd, modular inverse, and the Chinese Remainder Theorem — including the general non-coprime case with contradiction detection.

Background

Bézout identity: for any a, b there exist x, y with ax + by = gcd(a, b). Extended Euclid computes all three simultaneously: recursively from gcd(b, a mod b), or iteratively by carrying two coefficient rows through the ordinary gcd loop. O(log min(a, b)) steps.

Modular inverse: a⁻¹ mod m exists iff gcd(a, m) = 1; it’s the Bézout x of ax + my = 1. The Fermat shortcut a^(p−2) mod p works only for prime p; for composite m it’s ext_gcd or a^(φ(m)−1) — and φ needs m’s factorization, so ext_gcd is the default.

CRT: x ≡ a₁ (mod m₁), x ≡ a₂ (mod m₂) with coprime moduli has a unique solution mod m₁m₂. Non-coprime general case: let g = gcd(m₁, m₂); solvable iff g | (a₂ − a₁), and then unique mod lcm(m₁, m₂). k congruences: fold pairwise merges left to right.

Garner’s algorithm: reconstructs a big integer from its residues mod several primes in mixed-radix form, never materializing numbers larger than the answer — the standard finish for multi-modulus NTT (see phase-12 lab-07).

Interview Context

  • Codeforces / ICPC: modular inverse is used constantly (nCr mod p); CRT appears in mid-hard problems
  • Cryptography engineering: RSA keygen and RSA-CRT decryption are literally these primitives
  • Quant / research: exact rational arithmetic via multi-modular tricks
  • Systems: hashing, ring buffers, cycle alignment of periodic tasks
  • Standard FAANG: ext_gcd rare; pow(a, -1, m) occasionally in counting problems

When to Skip This Topic

Skip if any of these are true:

  • You don’t do modular-arithmetic problems at all (no competitive, no crypto)
  • You haven’t got plain gcd and fast exponentiation cold
  • Your only need is inverse mod a prime — memorize pow(a, p-2, p) and move on

Of the three Phase-13 number-theory labs, this one has the widest reach. Skip last.

Problem Statement

System of congruences, non-coprime moduli.

Given up to 10^5 congruences x ≡ aᵢ (mod mᵢ) with mᵢ ≤ 10^9 and NO coprimality guarantee, output the smallest non-negative x satisfying all of them, or “no solution”.

Constraints

  • 1 ≤ k ≤ 10^5 congruences, 1 ≤ mᵢ ≤ 10^9
  • Moduli arbitrary — shared factors and contradictions both possible
  • Combined modulus is lcm(m₁…mₖ): astronomically large — bigints required (Python native)
  • Wall-clock: < 2 sec

Clarifying Questions

  1. Are the moduli pairwise coprime? (No — the general case is the exercise.)
  2. Can the system be inconsistent? (Yes — detect and report “no solution”.)
  3. Are the aᵢ already reduced mod mᵢ? (Not guaranteed — reduce first; they may be negative.)
  4. How large does the answer get? (Up to lcm of all mᵢ — thousands of digits; fine in Python, needs bigint elsewhere.)

Examples

ext_gcd(240, 46) — full iterative table:

 r   240   46   10    6    4    2   0
 q          5    4    1    1    2
 s     1    0    1   −4    5   −9
 t     0    1   −5   21  −26   47
gcd = 2,  Bézout: 240×(−9) + 46×47 = −2160 + 2162 = 2

Three-congruence merge (Sunzi’s original):

x ≡ 2 (mod 3), x ≡ 3 (mod 5), x ≡ 2 (mod 7)
merge(2,3 | 3,5)  → x ≡ 8  (mod 15)
merge(8,15 | 2,7) → x ≡ 23 (mod 105)

Contradiction:

x ≡ 1 (mod 4), x ≡ 3 (mod 8): g = gcd(4,8) = 4, (3−1) % 4 = 2 ≠ 0 → no solution
(x ≡ 3 mod 8 forces x ≡ 3 mod 4)

Brute Force

Scan candidates x = 0, 1, …, lcm−1 and test all k congruences. lcm of 10^5 moduli up to 10^9 can exceed 10^100000. Dead on arrival.

Brute Force Complexity

  • Time: O(lcm × k) — unbounded for practical purposes
  • Space: O(1)

Optimization Path

  1. ext_gcd: gcd + Bézout coefficients in O(log min(a, b)).
  2. Two-congruence merge: x = a₁ + m₁t; need m₁t ≡ a₂−a₁ (mod m₂). With g = gcd(m₁, m₂): solvable iff g | (a₂−a₁); divide through by g, invert m₁/g mod m₂/g via Bézout, result unique mod lcm.
  3. Fold the merge over all k congruences; abort on first contradiction.
  4. Keep intermediates reduced (t mod m₂/g) so numbers grow no faster than the running lcm.

Final Expected Approach

def ext_gcd(a, b):
    old_r, r = a, b
    old_s, s = 1, 0
    old_t, t = 0, 1
    while r:
        q = old_r // r
        old_r, r = r, old_r - q * r
        old_s, s = s, old_s - q * s
        old_t, t = t, old_t - q * t
    return old_r, old_s, old_t

def modinv(a, m):
    g, x, _ = ext_gcd(a % m, m)
    if g != 1:
        raise ValueError("no inverse: gcd != 1")
    return x % m

def crt_merge(a1, m1, a2, m2):
    g, p, _ = ext_gcd(m1, m2)
    if (a2 - a1) % g:
        return None
    lcm = m1 // g * m2
    t = (a2 - a1) // g * p % (m2 // g)
    return (a1 + m1 * t) % lcm, lcm

def solve(congruences):
    a, m = 0, 1
    for ai, mi in congruences:
        merged = crt_merge(a, m, ai % mi, mi)
        if merged is None:
            return None
        a, m = merged
    return a

Starting from x ≡ 0 (mod 1) makes the fold uniform — every input congruence is just another merge.

Data Structures

  • Two integers (a, m): the running solution class x ≡ a (mod m)
  • No arrays, no recursion — the entire state is one residue class

Correctness Argument

  • ext_gcd invariant: old_s·a + old_t·b = old_r and s·a + t·b = r hold before and after every loop iteration (each row update is a linear combination of the previous two). At exit old_r = gcd with its Bézout pair.
  • Merge: x satisfies both congruences iff x = a₁ + m₁t and m₁t ≡ a₂−a₁ (mod m₂). This linear congruence in t is solvable iff g | (a₂−a₁); dividing by g gives (m₁/g)t ≡ (a₂−a₁)/g (mod m₂/g) with coprime coefficient, so t is unique mod m₂/g — hence x unique mod m₁·(m₂/g) = lcm. From m₁p + m₂q = g, p is exactly (m₁/g)⁻¹ mod (m₂/g).
  • Fold: by induction the running (a, m) encodes exactly the solution set of the congruences processed so far — either a full residue class mod lcm, or empty (early return). Final a ∈ [0, lcm) is the smallest non-negative solution.

Complexity

  • ext_gcd: O(log min(a, b)) ≈ ≤ 45 iterations for 10^9 inputs (worst case: Fibonacci pairs)
  • Full solve: k merges; bigint arithmetic on the running lcm of L bits makes it O(k · M(L)) — fast in practice because shared factors keep L from growing at every step
  • Space: O(L) for the running modulus

Implementation Requirements

  • Iterative ext_gcd (recursive dies at 10^5-scale inputs in Python only if misused, but iterative is the habit)
  • Reduce aᵢ mod mᵢ before merging — inputs may be negative or oversized
  • lcm as m1 // g * m2 — divide before multiplying
  • Reduce t mod m₂/g before forming x — keeps intermediates within ~2× the answer
  • Contradiction returns “no solution” immediately; don’t process remaining congruences
  • Final answer normalized to [0, lcm)

Tests

  • Sunzi: [(2,3), (3,5), (2,7)] → 23
  • Contradiction: [(1,4), (3,8)] → no solution
  • Redundant overlap: [(3,8), (3,4)] → 3 mod 8 (non-coprime, consistent)
  • Single congruence: [(5,7)] → 5; identity: [(0,1)] → 0
  • ext_gcd(240, 46) = (2, −9, 47); modinv(3, 7) = 5; modinv(4, 8) raises
  • Stress: random systems with lcm < 10^6 vs brute-force scan; assert identical including “no solution”
  • Performance: 10^5 random congruences, moduli ≤ 10^9, < 2 sec

Follow-up Questions

  • Why not Fermat a^(m−2) for the inverse? → That’s Fermat’s little theorem — prime modulus only. For composite m the exponent is φ(m)−1, and computing φ requires factoring m. ext_gcd needs neither primality nor factorization and is O(log m); it’s the default. (Python 3.8+: pow(a, -1, m) does this for you.)
  • This overflows in C++/Java — where? → m₁·m₂ up to 10^18 fits int64, but a₁ + m₁·t and (a₂−a₁)·p exceed it. Use __int128 (C++) or Math.multiplyHigh/mulmod (Java); for k-fold merges the running lcm needs a real bigint. Python is exempt — that’s why this lab is in Python.
  • What is Garner’s algorithm and when do you need it? → Given x mod p₁, …, x mod pₖ, Garner computes mixed-radix digits x = c₁ + c₂p₁ + c₃p₁p₂ + …, each cᵢ found with one modinv mod pᵢ. All arithmetic stays word-sized until final assembly — essential for reconstructing multi-modulus NTT results (phase-12 lab-07) in languages without cheap bigints.
  • Solve ax ≡ b (mod m) when gcd(a, m) = g > 1? → Solvable iff g | b. Then x ≡ (b/g)·(a/g)⁻¹ (mod m/g) — one solution mod m/g, hence exactly g distinct solutions mod m, spaced m/g apart.
  • Where does RSA use this? → Keygen: d = e⁻¹ mod φ(n) is one ext_gcd call. Decryption: RSA-CRT computes m mod p and m mod q with half-size exponents and moduli, then recombines with a precomputed q⁻¹ mod p — modmul is quadratic and exponents halve, so ~4× faster; every serious TLS stack ships it.
  • Non-coprime CRT vs coprime — what actually changes? → Coprime: g = 1, always solvable, modulus m₁m₂. Non-coprime: a solvability condition appears (g | a₂−a₁) and the combined modulus drops to lcm. The coprime formula silently gives wrong answers on non-coprime input — the classic buried bug.

Product Extension

  • TLS/SSH stacks: RSA-CRT private-key operations
  • Multi-modular exact linear algebra: solve mod several primes, CRT the result (SageMath, NTL)
  • Calendar/scheduling engines: alignment of periodic events is a congruence system
  • Distributed ID reconstruction: residue number systems for parallel arithmetic

Language/Runtime Follow-ups

  • Python: % returns non-negative, pow(a, -1, m) exists, bigints are native — the easiest language for this by a wide margin.
  • C++: % keeps the sign of the dividend — every reduce needs ((x % m) + m) % m; use __int128 for merge products.
  • Java: Math.floorMod, BigInteger.modInverse throws on non-coprime — matching our modinv contract.
  • Rust: i128 covers the pairwise merge; num-bigint for the fold.

Common Bugs

  1. Negative Bézout coefficient used raw: ext_gcd routinely returns negative x; every exit path needs a final % m normalization (which in C++ still isn’t enough — see sign bug).
  2. Coprime formula on non-coprime moduli: no crash, just wrong answers — the merge must go through g.
  3. Sign of (a₂−a₁) % g in C++: a negative nonzero remainder marks solvable systems as contradictions; Python’s % hides this class of bug entirely.
  4. lcm as m1 * m2 / g: overflows before the division in fixed-width languages; always m1 / g * m2.
  5. t not reduced mod m₂/g: correct value but intermediates balloon; in C++ this is an overflow, in Python a slowdown.
  6. aᵢ not pre-reduced: negative or oversized residues corrupt the solvability check.
  7. modinv without checking g: returns garbage “inverse” for gcd > 1; must raise or return sentinel.

Debugging Strategy

  • Assert g == ax + by after every ext_gcd call — one line, catches coefficient-row bugs instantly
  • Verify the final x against every input congruence before returning — O(k), cheap, self-verifying
  • Stress against brute-force scan on systems with lcm < 10^6, including engineered contradictions
  • Reproduce the 240/46 table by hand next to a printout of the loop state

Mastery Criteria

  • ext_gcd + modinv from memory in ≤ 10 min
  • General-case crt_merge + fold in ≤ 25 min
  • Reproduce the ext_gcd(240, 46) table by hand, no reference
  • State the non-coprime solvability condition and the merged modulus without hesitation
  • Explain why Fermat inverse fails for composite m in one sentence
  • Explain the RSA-CRT 4× speedup in under two minutes