Stern-Brocot Tree, in depth: mediants and the determinant invariant, paths as continued fractions, O(log N) oracle search, simplest fractions in an interval, and Farey sequences

Building the tree from mediants

Every node carries two bounds: its nearest ancestor to the left, a/b, and its nearest ancestor to the right, c/d. The node itself is the mediant (a+c)/(b+d). Its left child keeps the left bound and uses the node as the new right bound; its right child does the opposite. The root 1/1 is the mediant of 0/1 and 1/0.

The whole tree rests on one invariant: for the bounds of every node, bc − ad = 1. It holds for 0/1 and 1/0 (1·1 − 0·0 = 1), and it is preserved by each step, because replacing c/d with the mediant gives b(a+c) − a(b+d) = bc − ad = 1, and replacing a/b likewise. Three consequences follow:

• Ordering. bc − ad = 1 > 0 means a/b < c/d, and the mediant lies strictly between them. So the tree is a binary search tree by value.
• Lowest terms. The mediant m/n satisfies bm − an = 1 too, so any common divisor of m and n divides 1. Fractions in the tree never need reducing.
• Completeness. For any p/q in lowest terms, walking toward it from the root either hits it or narrows the interval; since a/b < p/q < c/d with bc − ad = 1 forces q ≥ b + d, the denominators grow and the walk must stop at p/q. Every positive rational appears exactly once.

Paths, matrices and Euclid

The walk is ordinary binary search with the mediant as the pivot. Comparing p/q with m/n by cross-multiplication keeps everything in integers.

from fractions import Fraction

def path_to(p, q):
    """L/R path from 1/1 to p/q (p, q > 0, coprime). Naive: O(depth)."""
    a, b, c, d = 0, 1, 1, 0            # left bound a/b, right bound c/d
    out = []
    while True:
        m, n = a + c, b + d            # the current node, the mediant
        if p * n == q * m:
            return "".join(out)
        if p * n < q * m:              # target lies left of the node
            c, d = m, n
            out.append("L")
        else:
            a, b = m, n
            out.append("R")

def run_length_path(p, q):
    """The same path straight from Euclid: R^a0 L^a1 R^a2 ..., last run minus one."""
    runs = []
    while q:
        runs.append(p // q)
        p, q = q, p % q
    runs[-1] -= 1
    return [("R" if i % 2 == 0 else "L", k) for i, k in enumerate(runs) if k]

L = ((1, 1), (0, 1))
R = ((1, 0), (1, 1))

def matmul(X, Y):
    return tuple(tuple(sum(X[i][k] * Y[k][j] for k in range(2)) for j in range(2))
                 for i in range(2))

def node_of(path):
    S = ((0, 1), (1, 0))               # columns: left bound 0/1, right bound 1/0
    for step in path:
        S = matmul(S, L if step == "L" else R)
    (n1, n2), (d1, d2) = S
    return Fraction(n1 + n2, d1 + d2)

assert path_to(3, 7) == "LLRR"
assert run_length_path(3, 7) == [("L", 2), ("R", 2)]
assert node_of("LLRR") == Fraction(3, 7)

Worked example, 3/7. Start between 0/1 and 1/0 at 1/1: 3/7 < 1/1, go L, bounds become 0/1 and 1/1, node 1/2. 3/7 < 1/2, go L: node 1/3. 3/7 > 1/3, go R: bounds 1/3 and 1/2, node 2/5. 3/7 > 2/5, go R: node 3/7. The path is LLRR.

Now compare with Euclid's algorithm on 3 and 7: quotients 0, 2, 3, so the continued fraction is [0; 2, 3]. The path is R⁰ L² R³⁻¹ = LLRR. In general the run lengths of the path are the continued-fraction terms, with one subtracted from the last. Continued fractions covers that expansion itself; here it tells us something about cost. The naive walk takes a₀ + a₁ + … steps, which is 1,000,000 steps for 1/1,000,000, while the run-length form takes only as many steps as Euclid's algorithm, O(log q).

The matrix form packs the two bounds into one 2×2 matrix, numerators on the top row and denominators on the bottom. Going left multiplies on the right by L, which keeps the left column and replaces the right column by the column sum; R does the reverse. Because the start matrix and both L and R have determinant ±1, the bc − ad = 1 invariant becomes 'the determinant never changes', and a run of k steps is a single multiplication by L^k = [[1, k], [0, 1]].

Searching with a comparison oracle

Suppose a value is a rational with denominator at most N, but all you can do is ask whether a given fraction is less than or equal to it. Typical cases are a hidden ratio inside a black-box system, a threshold found by a monotone test, or an interactive puzzle. Binary search on reals never terminates exactly; the Stern-Brocot tree does, and galloping along each run keeps it fast.

def gallop(ok):
    """Largest k >= 0 with ok(k) true; ok is monotone (true, then false) and ok(0) holds."""
    hi = 1
    while ok(hi):
        hi *= 2
    lo = hi // 2
    while hi - lo > 1:
        mid = (lo + hi) // 2
        if ok(mid):
            lo = mid
        else:
            hi = mid
    return lo

def find_hidden(le, N):
    """Recover a hidden positive rational x with denominator <= N.
    le(p, q) answers 'is p/q <= x?'. Invariant: a/b <= x < c/d."""
    a, b, c, d = 0, 1, 1, 0
    while b + d <= N:
        # run of R steps: largest k with (a + k c)/(b + k d) <= x
        k = gallop(lambda k: b + k * d <= N and le(a + k * c, b + k * d))
        a, b = a + k * c, b + k * d
        if b + d > N:
            break
        # run of L steps: largest k with (c + k a)/(d + k b) > x
        k = gallop(lambda k: d + k * b <= N and not le(c + k * a, d + k * b))
        c, d = c + k * a, d + k * b
    return Fraction(a, b)

Each run of identical moves is found with an exponential probe followed by binary search, costing O(log a_i) queries for a run of length a_i. The sum of log a_i is bounded by log q, and the number of runs is O(log q), so the whole search is O(log N) queries. In the test harness, 3,000 random targets with N = 10⁶, plus the adversarial cases 1/N, N/1 and the Fibonacci ratio 832040/514229, needed at most 66 oracle calls. The naive walk would need a million for 1/N. When the loop stops, a/b is the hidden value: the bounds are adjacent, so no fraction with denominator up to N lies strictly between them.

The simplest fraction in an interval

Because the tree is a search tree in which depth grows with denominator, the first node you meet that lies inside an interval is the fraction with the smallest denominator in that interval. This answers questions such as 'what is the simplest fraction that displays as 0.30 after rounding' or 'which small gear ratio lands inside a tolerance band'.

def smallest_between(x, y):
    """Fraction with the smallest denominator strictly inside (x, y), 0 <= x < y."""
    a, b, c, d = 0, 1, 1, 0
    while True:
        m = Fraction(a + c, b + d)
        if m <= x:      # jump right as far as possible while staying <= x
            k = gallop(lambda k: Fraction(a + k * c, b + k * d) <= x)
            a, b = a + k * c, b + k * d
        elif m >= y:    # jump left as far as possible while staying >= y
            k = gallop(lambda k: Fraction(c + k * a, d + k * b) >= y)
            c, d = c + k * a, d + k * b
        else:
            return m

assert smallest_between(Fraction(3, 10), Fraction(31, 100)) == Fraction(4, 13)
assert smallest_between(Fraction(314, 100), Fraction(315, 100)) == Fraction(22, 7)

The simplest fraction in (0.30, 0.31) is 4/13 ≈ 0.3077; the simplest in (3.14, 3.15) is 22/7. Both results were checked against a brute-force scan over denominators on 2,000 random intervals. Decide on open versus closed ends up front: with closed intervals, the endpoints themselves can be the answer, and the comparisons change from <= to <.

A close relative is the best approximation with a bounded denominator: walk toward x until the next node's denominator would exceed N, then return whichever bound is closer. Python ships this as Fraction.limit_denominator. CPython's implementation steps through continued-fraction convergents, which are exactly the run boundaries of the tree path, and then takes one last partial run (a semiconvergent) before comparing the two candidates. So Fraction(3.14159).limit_denominator(1000) returns 355/113 and Fraction(44100, 48000) reduces to 147/160: a 44.1 kHz to 48 kHz resampler interpolates by 160 and decimates by 147.

Farey sequences and the Calkin-Wilf tree

The Farey sequence F_n lists the fractions in [0, 1] with denominator at most n, in increasing order. It is exactly the in-order traversal of the Stern-Brocot subtree rooted at 1/2 with 0/1 and 1/1 added at the ends, pruned wherever the denominator would exceed n. The determinant invariant carries over: neighbours a/b, c/d in F_n satisfy bc − ad = 1, which the test harness confirmed for every n below 40.

def farey(n):
    out = [Fraction(0, 1)]
    def rec(a, b, c, d):
        if b + d > n:
            return
        rec(a, b, a + c, b + d)
        out.append(Fraction(a + c, b + d))
        rec(a + c, b + d, c, d)
    rec(0, 1, 1, 1)
    out.append(Fraction(1, 1))
    return out

assert [str(f) for f in farey(5)] == ["0", "1/5", "1/4", "1/3", "2/5", "1/2",
                                      "3/5", "2/3", "3/4", "4/5", "1"]

The Calkin-Wilf tree contains the same fractions on each level in a different order: each level is the bit-reversal permutation of the corresponding Stern-Brocot level. Its breadth-first order is generated by Stern's diatomic sequence, fusc(n)/fusc(n+1), which gives an O(1)-per-term enumeration of all positive rationals without repeats. Use Calkin-Wilf when you need to enumerate rationals, and Stern-Brocot when you need to search them by value.

Failure modes

• Unbounded naive walks. Walking one step at a time to 1/N or N/1 costs N iterations. If the input comes from users, that is a denial-of-service bug. Always take runs.
• Overflow in fixed-width integers. Cross-multiplying p · n against q · m doubles the bit width. In C, C++ or Rust use 128-bit products or check bounds before multiplying; in the galloping search, test b + k d ≤ N before computing the fraction so k cannot run away.
• Treating floats as exact. Fraction(0.1) is 3602879701896397/36028797018963968, the binary double, not 1/10. Convert decimal strings with Fraction("0.1"), or snap floats with limit_denominator and a sensible bound.
• Near-miss constants. 29.97 snapped with a bound of 2,000 gives 2997/100, but NTSC video is exactly 30000/1001 ≈ 29.97003. Snapping recovers the simplest nearby rational, not the true one; when a standard defines the ratio, use the standard's value.
• Signs, zero and infinity. The tree holds positive rationals only. Handle the sign separately, treat 0 as 0/1, and never let 1/0 escape as a result.

Choosing an approach

What to do next

1. Type in path_to and node_of and check that every fraction with p, q < 50 round-trips.
2. Implement find_hidden against a hidden Fraction and count oracle calls for 1/N, N/1 and a Fibonacci ratio.
3. Use smallest_between to choose display fractions for a set of rounded measurements.
4. Replace any float-based ratio detection in your code with limit_denominator and an explicit bound.
5. Read the Euclidean algorithm and the extended Euclidean algorithm to see why bc − ad = 1 is a Bézout identity.
6. Browse number theory for programmers for related integer tools.