Continued Fractions, in depth: Euclid's algorithm in disguise, convergents, best rational approximations, Pell equations and Wiener's attack

From division to an expansion

Take 415/93. Divide, keep the integer part, and flip the remainder:

415/93 = 4 + 43/93      -> a0 = 4, continue with 93/43
 93/43 = 2 +  7/43      -> a1 = 2, continue with 43/7
  43/7 = 6 +  1/7       -> a2 = 6, continue with 7/1
   7/1 = 7              -> a3 = 7, remainder 0, stop

415/93 = 4 + 1/(2 + 1/(6 + 1/7)) = [4; 2, 6, 7]

The pairs (415, 93), (93, 43), (43, 7), (7, 1) are exactly the steps of Euclid's algorithm on 415 and 93, and the integer parts 4, 2, 6, 7 are its quotients. So a rational number always has a finite continued fraction, its length is the number of Euclid steps, and that number grows only logarithmically in the denominator. An irrational number never terminates: the remainder is never zero, and for a real x the step is a = floor(x); x = 1 / (x - a).

One small ambiguity: a finite expansion can always end in two ways, because [..., a, 1] equals [..., a + 1]. The canonical form takes the last quotient greater than 1, and code that compares expansions must normalise first.

Convergents and the recurrence

Truncating the expansion after n terms gives the n-th convergent. You never evaluate the nested fraction; two running numerators and denominators suffice. With h(-1) = 1, h(-2) = 0, k(-1) = 0, k(-2) = 1:

h(n) = a(n) * h(n-1) + h(n-2)
k(n) = a(n) * k(n-1) + k(n-2)

The recurrence is the whole algorithm. It costs two multiply-adds per term, works left to right as quotients arrive, and never loses precision when the inputs are integers.

Code you can trust

Python's integers are arbitrary precision, which makes the exact versions short. The float version needs a stopping rule, because after roughly 15 to 17 significant digits the quotients describe rounding error, not the number.

from fractions import Fraction
import math

def cf_rational(p, q):
    """Partial quotients of p/q (q > 0). Same loop as Euclid's gcd."""
    out = []
    while q:
        a = p // q                 # floor division: correct for negative p too
        out.append(a)
        p, q = q, p - a * q
    return out

def convergents(quotients):
    """Yield (h, k) for each prefix of the expansion."""
    h0, k0, h1, k1 = 0, 1, 1, 0
    for a in quotients:
        h0, k0, h1, k1 = h1, k1, a * h1 + h0, a * k1 + k0
        yield h1, k1

def cf_real(x, max_terms=20):
    """Expansion of a float; use only the leading terms."""
    out = []
    for _ in range(max_terms):
        a = math.floor(x)
        out.append(a)
        frac = x - a
        if frac < 1e-12:
            break
        x = 1 / frac
    return out

assert cf_rational(415, 93) == [4, 2, 6, 7]
assert list(convergents([4, 2, 6, 7]))[-1] == (415, 93)

For a float input it is better still to convert it exactly with Fraction(x) and run cf_rational on the numerator and denominator. That expansion is exact for the binary value the float actually holds; you then decide how many leading terms are meaningful.

Why convergents are the best approximations

Convergents have four properties that make them useful, each provable from the recurrence by induction.

1. They are in lowest terms. h(n) k(n-1) - h(n-1) k(n) = (-1)^(n-1), so any common divisor of h(n) and k(n) divides 1.
2. They alternate. Even-indexed convergents sit below x, odd-indexed ones above, and each is closer than the one before.
3. The error is tiny. |x - h(n)/k(n)| < 1 / (k(n) * k(n+1)), and since k(n+1) >= k(n), it is below 1 / k(n)^2. A large next quotient means an unusually good convergent.
4. Legendre's converse. If a fraction p/q satisfies |x - p/q| < 1 / (2 q^2), then p/q is a convergent of x. This is the fact Wiener's attack uses.

Convergents are best approximations in a strong sense: no fraction with a smaller or equal denominator makes |q x - p| smaller. For the weaker question, the closest fraction with denominator at most L, the answer is either a convergent or a semiconvergent, a fraction of the form (m h(n) + h(n-1)) / (m k(n) + k(n-1)) with 0 < m < a(n+1). Python's Fraction.limit_denominator implements exactly this search.

π shows both. Its expansion begins [3; 7, 15, 1, 292, 1, 1, 1, ...], giving convergents 3, 22/7, 333/106, 355/113 and 103993/33102. The quotient 292 after 355/113 is why that fraction is famous: its error bound is divided by roughly 292, so it is accurate to about seven significant digits with a three-digit denominator. Asking Fraction(math.pi).limit_denominator(100) returns 311/99, which is not a convergent but the semiconvergent with m = 14 between 22/7 and 333/106: no convergent has a denominator between 7 and 106, and 311/99 beats 22/7. With a limit of 1000 the answer is 355/113.

Square roots and Pell equations

The square root of a non-square integer has a continued fraction that is eventually periodic, and the period can be computed with integers only, using the triple (m, d, a):

from math import isqrt

def sqrt_cf(n):
    """Return (a0, period) with sqrt(n) = [a0; period repeated]. n must not be a square."""
    a0 = isqrt(n)
    m, d, a, period = 0, 1, a0, []
    while a != 2 * a0:
        m = d * a - m
        d = (n - m * m) // d       # always exact
        a = (a0 + m) // d
        period.append(a)
    return a0, period

def pell(n):
    """Smallest positive solution of x^2 - n y^2 = 1."""
    a0, period = sqrt_cf(n)
    quotients = [a0] + period * 2  # two periods always suffice
    for x, y in convergents(quotients):
        if x * x - n * y * y == 1:
            return x, y

For n = 7, sqrt(7) = [2; 1, 1, 1, 4] with period 4, and the convergents reach 8/3, so 8^2 - 7 * 3^2 = 1. For n = 61 the period has 11 terms and the smallest solution is x = 1,766,319,049, y = 226,153,980; brute-force search over y would need hundreds of millions of steps, while the continued fraction needs about two dozen. Every other solution comes from powers of the smallest one, which is how Pell equations show up in counting problems and in some competitive programming tasks. Pell equations are Diophantine, like those in linear Diophantine equations, but quadratic, and the continued fraction is the standard efficient method for solving them.

Wiener's attack on small RSA exponents

RSA picks e d = 1 + k φ(N) for some integer k. Divide by d φ(N): e/φ(N) - k/d = 1/(d φ(N)). Since φ(N) is close to N, k/d is an extremely good approximation of the public ratio e/N. Wiener showed in 1990 that when d < N^(1/4) / 3 and the primes satisfy q < p < 2q, the approximation is good enough for Legendre's theorem, so k/d appears among the convergents of e/N. Each candidate is tested by recovering φ(N) and solving a quadratic for the primes:

from math import isqrt

def wiener(e, N):
    for k, d in convergents(cf_rational(e, N)):
        if k == 0 or (e * d - 1) % k:
            continue
        phi = (e * d - 1) // k
        s = N - phi + 1                      # p + q
        disc = s * s - 4 * N                 # (p - q)^2
        if disc >= 0 and isqrt(disc) ** 2 == disc:
            r = isqrt(disc)
            return d, (s + r) // 2, (s - r) // 2
    return None

# Toy key: p = 10007, q = 9973, d = 29 (deliberately tiny)
N, e = 99799811, 75695045
print(wiener(e, N))   # (29, 10007, 9973), found at convergent 22/29

The bound here is about 33, so d = 29 qualifies. The attack is why key generation should use a small public exponent such as 65537 and let d be whatever it comes out as, and never choose a small d to speed up decryption. CRT-based decryption, covered in RSA in depth, gets the speed safely.

Engineering ratios

Outside number theory, continued fractions answer "which small integer ratio is closest to this number?" A clock synthesiser that multiplies by M/D with both under 256, a gear train, a display aspect ratio, or a resampler converting between rates: each is a call to Fraction(target).limit_denominator(L), followed by a check that the numerator also fits. When the ratio is exactly rational the expansion terminates and returns it, for instance 44100/48000 reduces to 147/160.

Failure modes

• Trusting float tails. Quotients after the first dozen or so terms of a float expansion are rounding noise; a huge quotient late in the expansion usually marks where precision ran out.
• Truncating division with negatives. Languages whose integer division rounds toward zero produce wrong quotients for negative inputs. Use floor division.
• Overflow in the recurrence. Numerators and denominators grow at least like Fibonacci numbers; in 64-bit arithmetic check before multiplying or use big integers.
• Off-by-one seeds. Starting the recurrence with the wrong initial pairs shifts every convergent. Test against a known expansion such as 415/93.
• Ignoring the numerator bound. limit_denominator bounds only the denominator.
• Assuming convergents are enough. Under a denominator limit the best answer may be a semiconvergent, as 311/99 shows.

Trade-offs

Continued fractions are the fastest exact route to best approximations: logarithmic in the denominator, integer-only, with a proof of optimality. Walking the Stern-Brocot tree finds the same fractions one mediant at a time and is easier to explain, but takes as many steps as the sum of the quotients rather than their count, which is hopeless for π's 292. Binary search over Farey sequences is simple but slower still. For approximate engineering ratios, a brute-force scan over denominators is fine below a few thousand and lets you add arbitrary constraints; beyond that, use the expansion. The underlying modular toolkit is summarised in number theory for programmers.

What to do next

1. Implement cf_rational and convergents and check them against 415/93 = [4; 2, 6, 7].
2. Reproduce π's convergents and confirm that a limit of 100 gives 311/99 and a limit of 1000 gives 355/113.
3. Write sqrt_cf and pell; verify n = 7 gives (8, 3) and n = 61 the large solution.
4. Run the Wiener example, then generate a key with a large d and confirm the attack returns None.
5. Audit code that converts floats to ratios: use exact conversion and bound both numerator and denominator.
6. Add property tests: lowest terms, alternating sides, and the error bound for random rationals.