Rounding rational numbers using Farey/Cauchy sequence

A while back, for a project at work, I found I needed a way to round rational numbers. That is, given a number p/q, find the nearest number p'/q' such that both p' and q' are less than some limit N. I would be surprised if the algorithm I worked out is not already well known in some areas, but I couldn't find a description of it when I went looking, so I'm describing it here for the benefit of others.

Background

The algorithm makes use of some properties of rational numbers which can be observed from what is called a Farey sequence or from the Stern-Brocot tree. (It might be more accurate to attribute the sequence to one C. Haros, or to Cauchy, but J. Farey's name is the one that stuck.)

What makes Farey sequences useful here is that they provide a simple way to generate an ordered list of rational numbers in reduced form. Starting with the first sequence (consisting of 0/1 and either 1/1 or 1/0, depending on whether or not you think the latter is an abomination unto the rationals), you can generate the next Farey sequence by inserting between each pair of numbers its mediant. The tree structure implied by this is the Stern-Brocot tree, and so you can alternatively think of the Nth Farey sequence as the result of trimming the Stern-Brocot tree to a given depth and then performing an in-order traversal.

What is not immediately obvious is that the rational numbers generated in this way will be in numerical order and will be fully reduced. Farey noticed this property; Cauchy (and, earlier, Haros) proved it.

Algorithm

For convenience, I'll call the larger of p and q the size of the rational number p/q. (So, for example, 1/3 has size 3; 22/7 has size 22; etc.)

In general, I'm assuming that all numbers are in reduced form; that is, gcd(p,q)=1.

To restate the problem: You have a number, p/q, which you want to round to the nearest rational number p'/q' whose size is less than some limit N. You can do this by finding p/q on a Farey sequence, then searching forwards and backwards until you find the nearest adjacent numbers with a small enough size, and picking the "nearer" of the two.

Some subsidiary algorithms are useful.

Find-neighbors

Consider the first Farey sequence F(n) which contains a particular number p/q. The number has a left neighbor pL/qL and a right neighbor pR/qR. By the construction of the sequence, p = pL+pR and q = qL+qR, that is, p/q is the mediant of its neighbors and requires no reduction. The neighbors pL/qL and pR/qR can be easily calculated using the property of Farey sequences that qL p - pL q = 1 (and q pR - p qR = 1).

  1. Compute p1, the multiplicative inverse of p mod q; and q1, the inverse of q mod p. (These can be computed simultaneously by the extended Euclid's algorithm.)
  2. Since qL p = 1 (mod q), and -pL q = 1 (mod p), -pL and qL are equal (mod p and q) to the computed inverses.
  3. The mediant relation can be used to resolve the (mod p, q) ambiguity and to compute qR and pR.

The left and right numbers both have smaller or equal size to p/q. In addition, any later-numbered Farey sequences will be constructed from this one by inserting numbers whose size is at least as large as the numbers between which they are inserted. Therefore, a number's neighbors in the sequence in which it first appears are in fact the next and previous rational numbers of no greater size.

Rounding to the desired precision

The find-neighbors algorithm will find smaller-sized approximations to a given number, but those approximations may still be larger than the desired size. In many cases it will be necessary to apply find-neighbors more than once. It's also necessary to decide which of the two neighbors to use, either for further rounding or as the algorithm's result.

If find-neighbors produces two numbers of different size, both of which are larger than the desired size, then applying find-neighbors to the larger-sized result will produce the smaller result as one of its outputs. This can be seen both from the properties of find-neighbors, and from the fact that find-neighbors' results are adjacent members of the next-smaller Farey series. Repeated application of find-neighbors to the larger result will, at each step, shift either the right or left number farther from the initial number, without affecting the other one. Eventually, both the right and left results will be smaller than the desired size, and (again, by the definition of find-neighbors) will be the closest such rational numbers to the initial value.

However, if both the right and left numbers are larger than the desired size, find-neighbors can instead be applied to the smaller number; this will not alter the eventual result of the algorithm but will take fewer steps.

Once a pair of sufficiently small numbers is found, it's necessary to decide which one to use as the result of the rounding algorithm (that is, whether to round up or round down).

An obvious approach is to use the same rounding rules typically used for floating-point arithmetic, and to choose the result which is numerically closer to the original number.

Knuth (Art 2e 4.5.1) states that mediant rounding (choosing the right or left number according to whether the original number was above or below their mediant, instead of according to whether it was above or below their arithmetic mean) has been found to result in better numerical stability in long-running calculations.

It's also possible to avoid the question of rounding entirely and to retain both values, as a lower and upper bound on the computation's true value, and to carry this through future computations as in interval arithmetic.

Mediant rounding can be done particularly easily. Since the final pair of numbers will be the result of a single invocation of find-neighbors, and find-neighbors' input is the mediant of its outputs, the right or left number can be chosen based on whether the final find-neighbors was invoked on the left or right intermediate value. In fact, if at some point in the algorithm one number is small enough but the other isn't, the result of the final find-neighbors invocation doesn't have to be computed at all --- the result that is chosen will be the one that duplicates the pre-existing result from the previous iteration. As a result, if mediant rounding is desired, it's possible to use a simple "greedy" algorithm which terminates as soon as any acceptable result is computed:

  1. From the initial p/q, use find-neighbors to compute rational numbers L and R.
  2. As long as both L and R are larger than the desired size, apply find-neighbors to the smaller of the two in order to produce a new pair (L, R).
  3. If only one of (L, R) is small enough, terminate with that as the result.
  4. Otherwise, the result to choose depends on whether L or R was used as the most recent input to find-neighbors: if L was used, terminate with the resulting R; and vice versa.

For other forms of rounding, don't terminate early but perform the final step. Then choose L or R according to the desired rounding rule.

It should not be possible for L and R to be of equal size, except in the degenerate cases of 0/1, 1/1, and 1/0, since one of them will be the mediant of a pair including the other.

Complexity

My gut feeling is that the number of steps needed to round a number usually grows logarithmically with the amount of rounding. In the worst case, though, the algorithm is linear, reducing the size of the number only by 1 on each step for example when rounding integers (p/1) or fractions of the form 1/q.