René's URL Explorer Experiment

{"@context":"https://schema.org","@type":"DiscussionForumPosting","headline":"Document correctness proof for Fraction.limit_denominator","articleBody":"The correctness of the [`Fraction.limit_denominator`](https://github.com/python/cpython/blob/834064c19a110dad425dc290c91c0545eaa24471/Lib/fractions.py#L202) 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).\r\n\r\nFor 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.\r\n\r\nI'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.\r\n\r\nNote 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.\r\n\r\n## Set-up: simplifying the code\r\n\r\nWe'll prove correctness for a slightly simpler version of the code than what's currently in the `fractions` module. Here's that simpler version:\r\n\r\n```python\r\ndef limit_denominator(m: int, n: int, l: int) -\u003e tuple[int, int]:\r\n    \"\"\"\r\n    Given a fraction m/n and a positive integer l, return integers r and s such\r\n    that r/s is the closest fraction to m/n with denominator bounded by l.\r\n\r\n    m/n need not be in lowest terms, but n must be positive.\r\n\r\n    On return, 0 \u003c s ≤ l and gcd(r, s) = 1.\r\n    \"\"\"\r\n    a, b, p, q, r, s, v = n, m % n, 1, 0, m // n, 1, 1\r\n    while 0 \u003c b and q + a // b * s \u003c= l:\r\n        a, b, p, q, r, s, v = b, a % b, r, s, p + a // b * r, q + a // b * s, -v\r\n    t, u = p + (l - q) // s * r, q + (l - q) // s * s\r\n    return (r, s) if 2 * b * u \u003c= n else (t, u)\r\n```\r\n\r\nThe 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](https://github.com/python/cpython/blob/834064c19a110dad425dc290c91c0545eaa24471/Lib/fractions.py#L236-L237) 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.\r\n\r\n## Proof overview\r\n\r\nIn 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 line of code that defines $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 \u003c 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.\r\n\r\nIn more detail, we'll show in the sections below that after exiting the while loop and defining $t$ and $u$, the following hold:\r\n\r\n- $(ts - ru)v = 1$\r\n- $0 \u003c s \\le l$\r\n- $0 \u003c u \\le l$\r\n- $l \u003c s + u$\r\n- $r/s \\le m/n \u003c t/u$ if $v=1$, and $t/u \u003c m/n \\le r/s$ if $v = -1$\r\n\r\nWe now show that any fraction $y/z$ (with $z$ positive) strictly between $r/s$ and $m/n$ has $z \u003e l$. We have:\r\n\r\n$$z = 1z = (ts-ru)vz = (tz-yu)vs + (ys - rz)vu.$$\r\n\r\nNow $v$ is either $1$ or $-1$. If $v=1$ then $r/s \u003c t/u$, so $r/s \u003c y/z \u003c 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 \u003c r/s$, so $t/u \u003c y/z \u003c r/s$ and $ys - rz$ and $tz -yu$ are negative. So again $(ys - rz)v$ and $(tz - yu)v$ are positive integers.\r\n\r\nSo in either case, $(tz - yu)v \\ge 1$ and $(ys - rz)v \\ge 1$, hence\r\n$$z = [(tz-yu)v]s + [(ys - rz)v]u \\ge s + u \u003e l.$$\r\n\r\nSo 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.\r\n\r\nNote 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.\r\n\r\n## Details: loop invariants\r\n\r\nThe initial state, and the state after every iteration of the while loop, satisfies the following six invariants:\r\n\r\n- $(ps - rq)v = 1$\r\n- $ar + bp = m$\r\n- $as + bq = n$\r\n- $0 \\le b \u003c a$\r\n- $0 \\le q \\le s \\le l$\r\n- $0 \u003c s$\r\n\r\nThese 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$).\r\n\r\nIt also follows from these loop invariants that\r\n$$(pn - mq)v = a$$ and $$(ms - rn)v = b.$$\r\nTo 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$.\r\n\r\n## Details: loop termination\r\n\r\nFor 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$.\r\n\r\n## Details: after the loop\r\n\r\nWrite $k = \\lfloor (l - q) / s \\rfloor$, so that $t = p + kr$ and $u = q + ks$. We also define an integer $c = a - kb$. It's immediate from the definitions and the earlier loop invariants that:\r\n$$(ts - ru)v = 1$$\r\nand\r\n$$cr + bt = m$$\r\nand \r\n$$cs + bu = n$$\r\nand finally,\r\n$$(tn - mu)v = c.$$\r\n\r\nBy definition of floor, $k \\le (l-q)/s \u003c k+1$, so scaling by $s$ gives $q + ks \\le l \u003c q + ks + s$, hence\r\n$$u \\le l \u003c u + s.$$\r\n\r\nWe also note that since we've exited the loop at this point, either $b = 0$ or $0 \u003c b$ and $l \u003c q + \\lfloor a / b \\rfloor s$. In this second case, it follows that $\\lfloor (l - q) / s \\rfloor \u003c \\lfloor a / b \\rfloor$, hence that $k \u003c  \\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\r\n$$b \\le c.$$\r\nFurthermore, since $b$ is positive, it follows that\r\n$$0 \u003c c.$$\r\nMoreover, the above two inequalities also hold in the case $b=0$, since in that case $c = a - kb = a$ (and $0 \u003c a$, from our loop invariants). So $b \\le c$ and $0 \u003c c$ are always true after loop exit.\r\n\r\nWe'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 \u003c c = (tn - mu)v$: if $v = 1$ then these give $r/s \\le m/n \u003c t/u$, while if $v = -1$ then $t/u \u003c m/n \\le r/s$.\r\n\r\n## Tie-breaking, part I\r\n\r\nBy 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\r\n$$bu \\le cs.$$\r\nAdding $bu$ to both sides and using the fact that $cs + bu = n$, this is equivalent to the check\r\n$$2bu \\le n,$$\r\nwhich is the check that we use in the code.\r\n\r\nIn 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).\r\n\r\n## Tie breaking, part II\r\n\r\nThere'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.)\r\n\r\nFirst note that in general it's true that at any point before, during, or after the loop, we have:\r\n$$2bq \u003c n.$$\r\nSince $n = as + bq$, it's equivalent to show that $bq \u003c as$, but this follows from $bq \\le bs \u003c as$, since $s$ is positive and $b \u003c a$.\r\n\r\nNow in a tie-break situation we have $n = 2bu$, hence $2bq \u003c 2bu$, so $b$ must be positive and $q \u003c u$.\r\n\r\nNow 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 \u003c u$ and $q$ is nonnegative, it follows that $q = 0$.\r\n\r\nBut 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 \u003c u + s$.\r\n\r\n## Optimization: getting rid of the $0 \u003c b$ check\r\n\r\nSuppose that on entry to the `limit_denominator` function, we already know that $m/n$ is in lowest terms and that $l \u003c 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 \u003c n$.\r\n\r\nSo with $m/n$ already in lowest terms and $l \u003c 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 \u003c b$ check at the top of the while loop.\r\n","author":{"url":"https://github.com/mdickinson","@type":"Person","name":"mdickinson"},"datePublished":"2022-08-05T20:23:05.000Z","interactionStatistic":{"@type":"InteractionCounter","interactionType":"https://schema.org/CommentAction","userInteractionCount":8},"url":"https://github.com/95723/cpython/issues/95723"}