Extended Euclidean algorithm

The extended Euclidean algorithm runs the ordinary Euclidean algorithm for gcd(a, b) while carrying along the bookkeeping that answers a second question: which integer combination of a and b produces that gcd. It costs the same O(log min(a, b)) steps and it is the workhorse behind almost every "divide in modular arithmetic" move in this stack — most often disguised as a modular inverse, occasionally in its full form as the general solution of a linear equation in integers.

The algorithm

Euclid's algorithm repeatedly replaces (a, b) with (b, a mod b) until the second term is zero; the first term is then the gcd. The extended version threads two extra registers through the same loop, tracking how the current remainder is built out of the original a and b:

def egcd(a, b):
    old_r, r = a, b
    old_s, s = 1, 0
    while r:
        q = old_r // r
        old_r, r = r, old_r - q * r
        old_s, s = s, old_s - q * s
    return old_r, old_s          # g, and x with a*x ≡ g (mod b)

Every remainder in the sequence is an integer combination x·a + y·b, and the loop just applies the same subtraction to the coefficients that it applies to the remainders. When it terminates you have the gcd g together with a pair (x, y) satisfying a·x + b·y = g — the coefficients of Bézout's identity. The y register can be dropped when you only want x, since y is recoverable as (g - a·x) / b; that is why the version above tracks a single coefficient. Write it iteratively rather than recursively — the recursion is prettier on a blackboard but the loop has no stack cost and is the form both languages here use.

The common case: a modular inverse

Nine times out of ten the reason to reach for it is a modular multiplicative inverse: given a and m, find a⁻¹ with a·a⁻¹ ≡ 1 (mod m). Take b = m. If gcd(a, m) = 1, Bézout gives a·x + m·y = 1, and reading that modulo m collapses the second term to leave a·x ≡ 1 (mod m) — so x mod m is the inverse. Nothing extra is needed; the inverse is a by-product of the gcd computation.

Two consequences are worth internalising. First, the algorithm returns g as well as x, so it tells you whether the inverse exists rather than requiring you to check separately: g ≠ 1 means there is no inverse. Problem 0134 turns on exactly that boundary — it must end a fixed block of decimal digits while being divisible by the next prime, which needs 10ᵈ inverted modulo that prime, and the one excluded pair in the statement is precisely the one where the modulus shares a factor with 10ᵈ and the inverse does not exist. Second, an inverse is what lets you divide inside a modulus, which is why this technique shows up wherever a linear congruence has to be solved for its unknown, and as the glue step inside the Chinese Remainder Theorem. Problems 0160, 0932 and 0934 all use it that way: each merges congruences modulo coprime factors, and each merge is one inverse.

Problem 0271 shows the CRT glue in its other arrangement, and is worth seeing for that reason. It sums the cube roots of unity modulo a squarefree n with fourteen prime factors, and the roots modulo n are every combination of one root modulo each factor — hundreds of combinations, not one. Merging pairwise would repeat the inversions for every combination. Instead, each factor q gets one idempotent up front,

e_q = (n/q) · ((n/q)⁻¹ mod q)        ≡ 1 (mod q),  ≡ 0 (mod n/q)

after which any combination of per-factor residues r_q is just Σ r_q·e_q mod n — additions and multiplications only. The inverses are paid once per factor, not once per combination: a fixed basis, then cheap linear algebra. The same problem uses the inverse a second time, inside a Newton step: lifting a root of x³ − 1 from p to p^k by Hensel's lemma is ω ← ω − (ω³ − 1)·(3ω²)⁻¹, and the division by the derivative is, once again, an inverse modulo the current prime power — one that exists precisely because 3ω² is a unit when p ≠ 3, which is exactly why p = 3 is the case that needs separate handling.

There is a language asymmetry here that recurs in most of the problems below. Python has the inverse built in — pow(a, -1, m) since 3.8 — so the Python solution is a single expression and the technique is invisible in the source. The C port has no such primitive and carries a hand-written mod_inverse doing the loop above. When you compare a Python and a C solution in this stack and the only structural difference is a fifteen-line helper near the top of the C file, this is usually what it is. Problem 0250 is the asymmetry taken one step further: there the Python side needs no inverse at all, and only the C side does — see the next section.

When the divisor is not a unit

The g ≠ 1 verdict is not always the end of the road. Problem 0250 counts subsets by residue and needs whole rows of binomial coefficients modulo 10¹⁶, built by the multiplicative recurrence C(n, i+1) = C(n, i)·(n−i)/(i+1). That division is exact over the integers but not legal modulo 10¹⁶: half the denominators are even, a fifth are multiples of 5, and neither has an inverse. Python dodges the question by carrying the binomial as an exact big integer and reducing only when it stores it. The C port cannot, and the fix is worth knowing in general: split the running value into the part the modulus can see and the part it cannot,

C(n, i) = 2^e₂ · 5^e₅ · u        with gcd(u, 10) = 1

Strip every factor of 2 and 5 out of each numerator and denominator into the exponent counters e₂, e₅ (which may go up and down), and what remains of the denominator is coprime to 10¹⁶, so extended Euclid inverts it and the division becomes a multiplication on u. The coefficient is reassembled as u · 2^e₂ · 5^e₅ mod 10¹⁶ only when it is needed. The recipe generalises to any composite modulus: the primes it shares with your denominators are tracked as exponents, everything else is a unit, and units are exactly what the algorithm can invert.

Problem 0250 keeps its modulus and changes the number. Problem 0531 does the opposite: it keeps the number and changes the modulus. It is the Chinese Remainder Theorem with the coprimality taken away. For every pair n < m in a narrow window it solves x ≡ a (mod n), x ≡ b (mod m), and neighbouring moduli routinely share factors. Substituting x = a + n·t leaves n·t ≡ b − a (mod m), where n is not a unit. The linear congruence theorem covers this: with g = gcd(u, m), the congruence u·t ≡ c (mod m) is solvable exactly when g divides c, and dividing all three through by g leaves

(u/g)·t ≡ c/g   (mod m/g)        with gcd(u/g, m/g) = 1

So g ≠ 1 does not mean "no solution". It asks a sharper question, whether g divides c, and on a yes the inverse exists one level down, modulo m/g. In 0531 that one divisibility test decides whether a pair's system is consistent at all, and a consistent pair has a single solution class modulo lcm(n, m). The substitution also pays off at the end: with a in [0, n) and t in [0, m/g), the value a + n·t already lies in [0, lcm(n, m)), so it is the least solution with no further reduction. With some twelve million pairs, the per-pair cost matters too. Because gcd(n, m) = gcd(n, m − n) and m − n is below the window's width, each gcd starts from a small operand.

Division inside a closed form

Derivations often end in a fraction when the answer is wanted modulo something. The standard case is a geometric series far too long to sum term by term:

1 + r + r² + … + r^(n−1)  =  (r^n − 1) / (r − 1)

Problems 0335, 0479 and 0520 all reduce to this, at one modular power for the numerator and one inverse for the denominator. It is legitimate because the fraction is an integer: q·d = N survives reduction, so q ≡ N·d⁻¹ whenever d is a unit. Note what that condition does not ask for. The modulus need not be prime. Problem 0335 works modulo 7⁹, and its denominators are small numbers coprime to 7, which is all the inverse requires.

That is also the case for Euclid over the alternative. By Fermat's little theorem, a⁻¹ ≡ a^(p-2) (mod p), computable by exponentiation by squaring. It is a one-liner and worth knowing, but it needs a prime modulus and spends O(log p) multiplications where extended Euclid spends O(log m) cheap divisions. On 7⁹ it does not apply as written, and Euler's generalisation needs φ(m) in hand. Prefer Euclid unless you are already doing modular exponentiation anyway.

What the closed form does need is r − 1 to be a unit, so settle that before trusting it. There are three situations to tell apart:

  • The denominator is provably safe. This is the case to prefer, but the proof belongs in a comment. Problem 0479 has no branch because its denominator is a power of the summation index, which stays far below the prime modulus, so the inverse exists every time.
  • r ≡ 1 modulo a prime. The formula reads 0/0, while the sum is plainly the term count. Problem 0520 sums twenty-one series whose ratios run through small integers of both signs, so it needs a branch: ratio 1 returns the term count, reduced modulo p.
  • r − 1 shares a factor with a composite modulus. Take r = 8 modulo 7⁹: r ≢ 1, yet r − 1 has no inverse, so neither branch above helps. Widen instead: compute r^n modulo m·(r − 1), subtract 1, and divide exactly. None of these problems needs it, but it is the failure a prime modulus hides.

Problem 0520 also shows that the fraction need not wait until the last step. Its generating function expands cosh and sinh in eˣ, each bringing a ½, so the coefficients it multiplies are true fractions with denominators up to 1024, and only the final power sum is an integer. Rather than carry rationals, it fixes ½ as the residue 2⁻¹ mod p once at the start, and everything after is integer arithmetic. Compare 0250, where each intermediate was an integer and the trouble was the divisor. Here the intermediates are not integers at all, and that does not matter: reduction modulo p extends to every fraction whose denominator is a unit, so a rational simply is a residue. Any derivation over the rationals can run modulo p, provided p divides no denominator.

All of those divisions come out exact in the end. Problem 0731 needs the harder case, a floor. It reads a window of digits far into a Stoneham number, which splits term by term into quotients ⌊N/D⌋ mod M with N a power of ten too large to write down. Multiplying by D⁻¹ alone gives nonsense, because N/D is not an integer. The remainder is what repairs it:

N = q·D + r   ⟹   q ≡ (N − r)·D⁻¹ (mod M),     r = N mod D

This takes two moduli. Reduce N modulo D to find the remainder, then work modulo M with the remainder subtracted, which makes the division exact. Since N is a power, both reductions are single modular exponentiations. The identity needs only gcd(D, M) = 1, and D may be far larger than M, as it is in 0731. The remainders are not thrown away either: their fractions r/D add up to the small carry that the floor of the whole sum needs. Whenever you want the low digits of a quotient too large to form, this is the move.

The general case: a solution family

The inverse is a special case that discards most of what the algorithm computes. The full object is the general solution of a linear Diophantine equation c₁x + c₂y = R. It is solvable exactly when g = gcd(c₁, c₂) divides R; scaling the Bézout coefficients by R/g gives one particular solution (x₀, y₀), and every solution is

(x₀ + k·c₂/g,  y₀ − k·c₁/g)     for k ∈ ℤ

— a one-parameter family, a lattice line. That shift vector is the useful part: it converts "does an integer solution exist in this box?" into two linear inequalities on k, one from each variable's interval, and the answer is whether the two k-ranges intersect. Problem 0236 (the Simpson's paradox hamper table) leans on this directly: its feasibility test reduces to a three-variable equation with each variable confined to a known interval, and after iterating the narrowest interval the remaining two-variable equation is decided in O(1) per candidate by exactly this range-intersection test. A brute-force scan over the widest interval instead would have been thousands of times the work.

Between the two cases sits a pattern common enough to name: count the solutions of a linear congruence in a range rather than find one. The solutions of a·x ≡ c (mod m) with gcd(a, m) = 1 are the single class x ≡ c·a⁻¹ (mod m), so the number of them in [lo, hi] is an arithmetic-progression count — O(1) once you have the inverse, however wide the range. Problem 0198's inner count has exactly this shape, with r an inverse that falls out of the Farey-neighbour relation bc − ad = 1 — a Bézout identity the problem hands you rather than one you have to compute. Problem 0251 (Cardano triplets) shows it at scale: after splitting b against n by their gcd, the divisibility condition becomes 8·n₁·d ≡ 3 (mod b₁²), so for each coprime pair (n₁, b₁) the admissible d are one residue class modulo b₁², its smallest member is 3·(8n₁)⁻¹ mod b₁², and the count up to a bound is a single floor division. Note how existence is settled by structure rather than tested: b₁ is forced odd, and pairs with gcd(n₁, b₁) > 1 are skipped before the inverse is ever asked for. A search over triplets collapses into one inverse per (n₁, b₁) pair.

Getting it right

  • Normalise the sign. The Bézout x can be negative; the inverse must land in [0, m). Python's % already returns a non-negative result for a positive modulus, so x % m suffices. C's % follows the sign of the dividend, so the idiom is ((x % m) + m) % m — and the same correction is needed on any other quantity you reduce, not just the inverse.
  • Watch the coefficient growth. The remainders shrink monotonically but the coefficients grow; they stay bounded by roughly max(a, b), so the danger is not the coefficients themselves but the x·a products you form afterwards. Reduce operands modulo m before multiplying, and reach for a wider type when the intermediate does not fit — problems 0236, 0271, 0731 and 0932 all do the arithmetic in __int128 on the C side for this reason, and 0250, whose modulus 10¹⁶ already fills most of a 64-bit word, runs the Euclid loop itself in __int128, where Python's arbitrary-precision integers simply absorb it.
  • Floor, don't truncate. The recurrence and the interval clipping both assume floored division. C truncates toward zero, so a C port that clips a solution family to an interval needs an explicit floored-division helper; getting this wrong shows up as an off-by-one only on negative inputs, which is exactly the case a quick test is least likely to cover.
  • Recognise the shape. If a derivation ends in "…so k ≡ something (mod n)" or "…so this must be a multiple of M", you are one inverse away from a closed form and should stop searching. That is the real payoff: the technique routinely replaces a loop over candidates with a single expression.

Problems14 / 14

  • ● 0134 — Prime Pair Connection
  • ● 0160 — Factorial Trailing Digits
  • ● 0198 — Ambiguous Numbers
  • ● 0236 — Luxury Hampers
  • ● 0250 — $250250$
  • ● 0251 — Cardano Triplets
  • ● 0271 — Modular Cubes, Part 1
  • ● 0335 — Gathering the Beans
  • ● 0479 — Roots on the Rise
  • ● 0520 — Simbers
  • ● 0531 — Chinese Leftovers
  • ● 0731 — A Stoneham Number
  • ● 0932 — $2025$
  • ● 0934 — Unlucky Primes