The dedoublifier, part four

I said last time that binary-searching the rationals (WOLOG between zero and one) for a particular fraction that is very close to a given double does not really work, because we end up with only fractions that have powers of two in the denominators. We already know what fraction with a power of two in the denominator is closest to our given double; itself!

But there is a way to binary-search the rationals between zero and one, as it turns out. Before we get into it, a quick refresher on how binary search works:

The idea is that we have a target, and we are in a loop searching for a match to the target. We have upper and lower bounds. Each time through the loop we find an element that is between the upper and lower bounds; if it is a match to the target, we're done. If not, then the new element becomes either the new upper bound or the new lower bound. In this way every time through the loop we have a more and more precise bound, and so we get closer and closer to the matching element.

The tricky bit is the part that I did not specify in my algorithm sketch just now: given an upper and a lower bound, how do you find the element "in the middle"? We tried taking the mean of the upper and lower bounds but saw that it doesn't work.

The solution for our problem is to not use the mean but rather the mediant. If we have two non-negative fractions in simplest form, a/c and b/d then their mediant is (a+b)/(c+d). This is also known colloquially as the "freshman sum" because young math students often forget that this is not how you add fractions. We'll also need their determinant which is the integer bc - ad.

You might not have heard of the mediant before; it has many interesting properties that are germane to our interests. (I must emphasize again that the mediant is only well-defined on non-negative fractions in their simplest form; throughout let's assume that all fractions I'm considering are not only in their simplest form, but also between zero and one.)

Let's suppose we have fractions q = a/c and r = b/d, and their mediant m = (a+b)/(c+d). Let's also assume that q is less than r, and both are between 0 and 1 inclusive.

• m is strictly between q and r. The easy proof is left as an exercise.
• If the determinant of q and r is 1, then the determinant of q and m is 1, and the determinant of m and r is also 1. Again, this is an easy proof.
• If q and r have determinant 1 then m is the fraction between q and r that has the smallest denominator. For example, suppose we have 1/4 and 1/3. Their mediant is 2/7. There is no fraction between 1/4 and 1/3 with a smaller denominator in simplest form. 1/6 and 1/5 are too small and 2/5 and 5/6 are too large. This is a little harder to prove; see the Wikipedia page for a sketch of the proof and a proof of the converse.

Now it should be clear why the mediant is the right choice when binary-searching the rationals in an attempt to discover a simple fraction:

• It divides and conquers. Just like the mean, the mediant is always between the upper and lower bounds.
• If we start at 0/1 and 1/1 as the bounds then the mediant is always the fraction between the bounds with the lowest denominator. (Do you see why? Hint: the determinant of 0/1 and 1/1 is 1.)

Let me make sure the implications of this are clear. Using the mediant as our dividing point does not produce a sequence of always-better approximations. (Neither, for that matter, does using the mean; both produce sequences of better and better bounds, but the midpoint chosen need not always be closer to the target than the previous midpoint.) Rather, it produces a sequence of the best approximations that have small denominators. The denominators, being powers of two, get extremely large extremely quickly when using the mean. When using the mediant, we are guaranteed that the sequence of approximations has the smallest possible denominators.

I think at this point it's best if we write some code.

public static IEnumerable<Rational> Approximations(this Rational x){ Rational fraction = x.AbsoluteValue.FractionalPart; if (fraction == 0) { yield return x; yield break; } Rational whole = x.AbsoluteValue.WholePart; Debug.Assert(fraction > 0); Debug.Assert(fraction < 1); Rational lower = 0; Rational upper = 1; while(true) { Rational mediant = Rational.FromIntegers( lower.Numerator + upper.Numerator, lower.Denominator + upper.Denominator); yield return x.Sign * (whole + mediant); if (fraction < mediant) upper = mediant; else if (fraction > mediant) lower = mediant; else break; }}

The code should be pretty straightforward. We simplify the problem by only considering fractions strictly between 0 and 1; if that's not what we were given then we extract the fractional part from the input and work with that. Every time through the loop we yield the mediant of our two bounds; if we happened to exactly hit the rational given then we cannot make any better approximation.

Let's take look at it in action.

Console.WriteLine(string.Join(", ", Rational.FromDouble(5.0/7.0).Approximations().Take(10)));

This gives us:

1/2, 2/3, 3/4, 5/7, 8/11, 13/18, 18/25, 23/32, 28/39, 33/46

Now the truth of what I said earlier is perhaps easier to see. This is not a sequence of ever-better approximations; 5/7 is by far the best approximation on this list and each element on the rest of this list will be worse for a long, long time. The denominator must get up into the quadrillions before we start getting better approximations, remember. Rather, the approximations after it are the best approximations that have slightly larger denominators.

How then can we get the desired target, 5/7, out of here? The sequence operators let us represent the desired semantics quite easily. We could use TakeWhile to say "stop when the approximation is within some epsilon of the target". Or we could use OrderBy to sort out the best of a set of approximations:

var r = Rational.FromDouble(5.0 / 7.0);Console.WriteLine(string.Join(", ", r.Approximations() .TakeWhile(x => x.Denominator <= 1000) .OrderBy(x => (x - r).AbsoluteValue) .Take(5)));

Which produces

5/7, 713/998, 708/991, 703/984, 698/977

There are any number of other such heuristics you could use.

This technique could also be used to find low-denominator approximations of irrational numbers. What is the fraction with denominator under 1000 that is the best approximation for pi?

var r = Rational.FromDouble(Math.PI);Console.WriteLine( r.Approximations() .TakeWhile(x => x.Denominator <= 1000) .OrderBy(x => (x - r).AbsoluteValue) .First());

Turns out that 355/113 is way better than 22/7. Who knew?

For more on binary-searching the rationals, see the Wikipedia page on the Stern-Brocot tree, which is the "binary search tree" implied by using the mediant. There is also some interesting material in there on the continued-fraction interpretation of the tree, which I decided to not go into in this series.

Well that's it for FAIC for 2015. Have a festive holiday season everyone and we'll see you for more fabulous adventures in 2016!