Document correctness proof for Fraction.limit_denominator #95723

The correctness of the Fraction.limit_denominator method is not immediately obvious, and extracting a proof of correctness from external sources is nontrivial: Wikipedia has statements without proofs at https://en.wikipedia.org/wiki/Continued_fraction#Best_rational_approximations (and at the time of writing those statements aren't 100% correct), and a typical elementary number theory book will have you wade through a whole chapter of theory on continued fractions to get to something useful. In spite of that, a targeted correctness proof for limit_denominator is straightforward, and doesn't require developing the whole theory of continued fractions (or even mentioning them).

For future maintainers of Fraction.limit_denominator, it would be useful to have the correctness proof recorded somewhere: when I recently had to update the limit_denominator code in #93730, I had to stare at the code for a good hour to figure out why it worked at all, despite the fact that I wrote that code in the first place (albeit a good number of years ago). Moreover, the new tie-break condition, while more efficient than the old one, is also less obvious.

I'm not sure what the best place to record the proof is: the fractions.py source? A text file alongside the source? For now, I'll record the proof in this issue, which is at least better than nothing.

Note on originality: the proof below is my own, and not a transcription from another source. I haven't seen any proof along these lines in the literature, but this is an extremely well-traveled road, so it's very likely that there's something closely resembling it somewhere. References would be welcome.

Set-up: simplifying the code

We'll prove correctness for a slightly simpler version of the code than what's currently in the fractions module. Here's that simpler version:

def limit_denominator(m: int, n: int, l: int) -> tuple[int, int]:
    """
    Given a fraction m/n and a positive integer l, return integers r and s such
    that r/s is the closest fraction to m/n with denominator bounded by l.

    m/n need not be in lowest terms, but n must be positive.

    On return, 0 < s ≤ l and gcd(r, s) = 1.
    """
    if l < 1:
        raise ValueError("max_denominator should be at least 1")
    
    a, b, p, q, r, s, v = n, m % n, 1, 0, m // n, 1, 1
    while 0 < b and q + a // b * s <= l:
        a, b, p, q, r, s, v = b, a % b, r, s, p + a // b * r, q + a // b * s, -v
    k = (l - q) // s
    t, u = p + k * r, q + k * s
    return (r, s) if 2 * b * u <= n else (t, u)

The differences with the existing code are mostly cosmetic: there are some de-optimizations (like computing a//b several times); we deal purely with integers to avoid the Fraction unpacking and repacking; some variables are renamed, and we've removed the fast path where a Fraction is returned unaltered if its denominator is already smaller than l. The most significant difference is that we've unrolled the first loop-and-a-half of the while loop, with the side-effect of shifting the loop exit condition to the top of the loop. There's also an extra variable v, which alternates between values 1 and -1. It's not needed for computing the final result, but we'll be making use of it in the statements and proofs below. Finally, the original version of the code doesn't check that b is positive before entering the loop - that check proves to be unnecessary if the fast path is present.

Proof overview

In a nutshell, the main loop is a variant of the extended Euclidean algorithm applied to $m$ and $n$, with an extra termination condition designed to make sure that the denominators don't go out of bounds. Together with the following lines of code that define $t$ and $u$, the main loop produces fractions $r/s$ and $t/u$ which bracket the target fraction $m/n$, and whose denominators are small (i.e., no greater than the limit $l$). With just two additional pieces of information, namely that $ts - ru = ±1$ and that $l < s + u$, we can deduce that all fractions strictly between $r/s$ and $t/u$ have denominator greater than $l$, hence that either $r/s$ or $t/u$ is the fraction that we want. The final line of the code then determines which of those two fractions is closest to $m/n$ and returns it, returning $r/s$ in the case of a tie-break.

In more detail, we'll show in the sections below that after exiting the while loop and defining $t$ and $u$, the following hold:

• $(ts - ru)v = 1$
• $0 < s \le l$
• $0 < u \le l$
• $l < s + u$
• $r/s \le m/n < t/u$ if $v=1$, and $t/u < m/n \le r/s$ if $v = -1$

We now show that any fraction $y/z$ (with $z$ positive) strictly between $r/s$ and $t/u$ has $z > l$. We have:

$$z = 1z = (ts-ru)vz = (tz-yu)vs + (ys - rz)vu.$$

Now $v$ is either $1$ or $-1$. If $v=1$ then $r/s < t/u$, so $r/s < y/z < t/u$. So $ys - rz$ and $tz - yu$ are positive, hence $(ys - rz)v$ and $(tz - yu)v$ are positive integers. Conversely if $v=-1$ then $t/u < r/s$, so $t/u < y/z < r/s$ and $ys - rz$ and $tz -yu$ are negative. So again $(ys - rz)v$ and $(tz - yu)v$ are positive integers.

So in either case, $(tz - yu)v \ge 1$ and $(ys - rz)v \ge 1$, hence
$$z = [(tz-yu)v]s + [(ys - rz)v]u \ge s + u > l.$$

So either $r/s$ or $t/u$ must be the closest fraction to $m/n$ amongst all fractions with small (i.e., bounded by $l$) denominator, and we merely have to determine which one is closer. The details of this are given in the "Tie-breaking" section below. The proofs of the five properties that we used above are in the "Details" sections below.

Note that it follows from $(ts - ru)v = 1$ that both $\gcd(t, u) = 1$ and $\gcd(r, s) = 1$, so whichever fraction is returned is already normalised. The original limit_denominator code makes use of this to avoid doing extra unnecessary gcd computations when creating the returned Fraction instance.

Details: loop invariants

The initial state, and the state after every iteration of the while loop, satisfies the following six invariants:

• $(ps - rq)v = 1$
• $ar + bp = m$
• $as + bq = n$
• $0 \le b < a$
• $0 \le q$
• $0 < s \le l$

These are all easily proved by induction: they hold for the initial state by inspection, and it's easy to see that the way that the variables are updated within the while loop keeps all of these properties invariant. Note that at the end of an iteration of the while loop, the condition $s \le l$ is exactly the condition that the loop was entered (that $q + \lfloor a / b \rfloor s \le l$).

It also follows from these loop invariants that
$$(pn - mq)v = a$$ and $$(ms - rn)v = b.$$
To see these, simply expand the values of $m$ and $n$ using the second and third invariants, and then collapse using the first. For example: $(pn - mq)v = p(as + bq)v - (ar + bp)qv = (ps - rq)va = a$.

Details: loop termination

For a complete proof of correctness we need to establish that our while loop is not infinite - that it terminates after a finite number of iterations. But this follows immediately from the observation that the value of $b$ is nonnegative and that it strictly decreases at every iteration (because the new value of $b$ at each iteration is computed from the old values as $a \bmod b$). Thus the total number of iterations can never be more than the initial value $m \bmod n$ of $b$.

Details: after the loop

We have $t = p + kr$ and $u = q + ks$, with $k = \lfloor (l - q) / s \rfloor$. We also define an integer $c = a - kb$. It's immediate from the definitions and the earlier loop invariants that:
$$(ts - ru)v = 1$$
and
$$cr + bt = m$$
and
$$cs + bu = n$$
and finally,
$$(tn - mu)v = c.$$

By definition of floor, $k \le (l-q)/s < k+1$, so scaling by $s$ gives $q + ks \le l < q + ks + s$, hence
$$u \le l < u + s.$$
Since also $s \le l$ (from the loop invariants) we have $s \le l < u + s$, from which $0 < u$.

We also note that since we've exited the loop at this point, either $b = 0$ or $0 < b$ and $l < q + \lfloor a / b \rfloor s$. In this second case, it follows that $\lfloor (l - q) / s \rfloor < \lfloor a / b \rfloor$, hence that $k < \lfloor a / b \rfloor$, or equivalently $k+1 \le \lfloor a / b \rfloor$. So $(k + 1)b \le a$ and $b \le a - kb$. So
$$b \le c.$$
Furthermore, since $b$ is positive, it follows that
$$0 < c.$$
Moreover, the above two inequalities also hold in the case $b=0$, since in that case $c = a - kb = a$ (and $0 < a$, from our loop invariants). So $b \le c$ and $0 < c$ are always true after loop exit.

We've now established almost all the facts that were used in the proof overview above. The only thing that we're missing is that $m/n$ lies between $r/s$ and $t/u$ (and can be equal to $r/s$, but not to $t/u$). But this follows from the facts that $0 \le b = (ms - rn)v$ and that $0 < c = (tn - mu)v$: if $v = 1$ then these give $r/s \le m/n < t/u$, while if $v = -1$ then $t/u < m/n \le r/s$.

Tie-breaking, part I

By the time we get to the final line of limit_denominator, we know that either $r/s$ or $t/u$ is our best approximation; we need to choose the one that's closest to $m/n$. So we want to compare $|m/n - r/s|$ with $|t/u - m/n|$. Scaling up by $nsu$, we need to compare $|(ms - rn)u|$ with $|(tn - mu)s|$. Since $v$ is either plus or minus one, we can insert a factor of $v$ without affecting the absolute value, so we're now comparing $|(ms-rn)vu|$ with $|(tn - mu)vs|$. But from above, $(ms-rn)v = b$ is nonnegative, as is $(tn - mu)v = c$. So the comparison we want to make is
$$bu \le cs.$$
Adding $bu$ to both sides and using the fact that $cs + bu = n$, this is equivalent to the check
$$2bu \le n,$$
which is the check that we use in the code.

In the case that $r/s$ and $t/u$ are equidistant from $m/n$, we have $bu = cs$ (or equivalently $2bu = n$) and we end up returning the bound $r/s$. The comments in the original limit_denominator implementation claim that the bound with smaller denominator is returned in that case: i.e., that $s \le u$. But $cs = bu \le cu$ (since $u$ is positive and $b \le c$), so $s \le u$ (since $c$ is positive).

Tie-breaking, part II

There's a second part to the tie-breaking comment in the limit_denominator source, that in the case of a tie, if both denominators are equal, then the lower bound is returned. (This can only happen in extremely limited circumstances, namely when $l = 1$ and $m/n$ is exactly halfway between two integers.)

First note that in general it's true that at any point before, during, or after the loop, we have:
$$2bq < n.$$
Since $n = as + bq$, it's equivalent to show that $bq < as$, but this follows from $bq \le bs < as$, since $s$ is positive and $b < a$.

Now in a tie-break situation we have $n = 2bu$, hence $2bq < 2bu$, so $b$ must be positive and $q < u$.

Now suppose that the denominators of our two candidates are equal: $s = u$. Since $(ts - ru)v = 1$, $s$ and $u$ are relatively prime, and since they're also positive, the only possibility is that $s = u = 1$. Now since $q < u$ and $q$ is nonnegative, it follows that $q = 0$.

But an examination of the while loop shows that the only time $q = 0$ is before entry to the while loop: after one or more iterations of the while loop, $q$ is always positive. So the only way to reach this situation is if the body of the while loop must have executed zero times, and $v = 1$. Now substituting back into $(ts - ru)v = 1$ gives $t = r + 1$, so $r/s = r$ is the lower bound for $m/n$ and $t/u = t$ is the upper bound; $m/n$ lies exactly halfway between two consecutive integers. Note also that in this case we must have $l = 1$, since $l < u + s$.

Optimization: getting rid of the $0 < b$ check

Suppose that on entry to the limit_denominator function, we already know that $m/n$ is in lowest terms and that $l < n$. Then on exit of the while loop, $b$ must be positive. For if $b = 0$ then the loop invariant $ar + bp = m$ gives $ar = m$, while the loop invariant $as + bq = n$ gives $as = n$. So since $m$ and $n$ are relatively prime, $a = 1$, $m = r$ and $n = s$. But $s \le l$ (another loop invariant from earlier), so $n \le l$, a contradiction to our assumption that $l < n$.

So with $m/n$ already in lowest terms and $l < n$, $b$ is positive at loop exit time, and since it only decreases during the loop, it's positive all the way through the calculation. So in that case, there's no need to keep the $0 < b$ check at the top of the while loop.