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!

Advertisements

19 thoughts on “The dedoublifier, part four

  1. Interesting series! I once had to implement from sratch a Rational Number implementation with all kind of yucky features like supporting NaN, Infinity, etc. Among other stuff, I had to provide a clever way of giving a rational approximation of any given real number. The path I took was using Continued Fractions, Allthough a whole lot clunkier, it worked out pretty well for our needs. If anyone is curious, you can find it at RationalNumber

  2. The mediant was more interesting the problem itself.

    BTW, I would try to solve this by searching for repeating bits, but then I realized that numbers like 1/67 starts repeating only after 67 binary digits and doubles only have 52.

  3. The Haskell implementation of approxRational is also interesting I think:

    https://hackage.haskell.org/package/base-4.8.1.0/docs/src/Data.Ratio.html#approxRational

    Internally, it effectively constructs a continued fraction approximation to a given number. This has the advantage of converging very quickly: each approximation is guaranteed to be better than the previous one. For example, for 5/7 (as a double), approxRational outputs 5/7, and then 6433713753386423 / 2^53.

    • Yeah, I considered going with a continued fraction approach, but really I wanted to talk about the binary searching algorithm here. Maybe I’ll do a bit on the continued fraction approach next year.

  4. Pingback: The Morning Brew - Chris Alcock » The Morning Brew #1988

  5. Thinking about your Approximation-of-pi function, you didn’t write it that way but I could imagine that enumerable object going into an indefinitely producing state. Its fine if you use Take/TakeWhile (as your example did) but naively call ToList and you’ve got a weird out-of-memory error.

    Maybe if there was a way for the caller of an enumerable to indicate “I have no terminating condition except your MoveNext returning false.” to which the enumerable would respond “My MoveNext will always return true, throw exception.” But for that to work, the indication would need to passed along by Select/Where/etc and some compiler cleverness to figure out if a foreach has a break/return statement inside. IE, probably not worth it just to avoid an easily diagnosed error.

      • This would not require a halting-problem solution.

        My idea was for the enumerable to receive some indication that it was being looped in a way that will probably never stop, such as a foreach loop with no break/return/goto statements or a .ToList call. The author of the enumerable probably knows if it will loop indefinitely or not so can make a check that the code on the other end has a separate terminating condition.

        Its a bad idea for so many reasons. I raise it because I’m interested in the implications of an enumerable with a MoveNext that never returns false.

    • Yeah, I think this is something we just have to live with when working with IEnumerable. It’s designed to allow effectively infinite sequences. When using it, you surely have some idea of what it represents and whether or not that sequence could be infinite. If it could (as in the case of approximating an irrational number), you just have to be ready to stop it manually.

  6. Pingback: Dew Drop – December 11, 2015 (#2150) | Morning Dew

  7. In the interest of simplicitly, one IEnumerable design is overloaded for two different purposes:

    1. There is a collection of items which may be read out sequentially.

    2. There is a means of sequentially generating and returning items as many items as a caller wants, at least up to some limit.

    The lines between the two purposes aren’t always clear, since there are some objects of the second type which can generate a finite number of items faster than objects of the first type would be able to read them out, but there are some kinds of consumer algorithms that are only suitable for use with implementations of the first type and some consumer algorithms that are better suited for implementations of the second (I don’t think that there are any consumer algorithms that are especially unsuitable for implementations of the first type, but there are some tasks which algorithms restricted to the first type can perform on the first type more efficiently than could algorithms that could work with both types).

    It’s too bad there’s no standard way for an implementation of IEnumerable to give consumers any clue as to how it should best be treated. For example, some implementations of IEnumerable guarantee that every future enumeration can promise that all future enumerations will always yield the same sequence, regardless of who has a reference to the object being enumerated; some can guarantee that the thing being enumerated will always yield the same sequence unless someone with a reference to it modifies it; some may return different data on every enumeration. An algorithm which needs to be able to iterate something multiple times may need to copy it to an array if the item can’t guarantee that repeated iterations will yield the same data, but such copying may be wasteful if repeated iterations would yield the same data.

  8. It looks like the algorithm is only converging linearly in the 5/7 example. Because 5/7 is such a good approximation, it is continually chosen as the lower bound. As a result the new denominator is only increasing by 7 each time. While you could argue it has already converged very closely to the answer, it will still take O(N) steps to reach the exact value.

    This is a point in favour of continued fractions / euclidean algorithm. The worst case is that all coefficients are 1, which still means the denominator will increase exponentially in the golden ratio.

  9. Pingback: 주간닷넷 2015년 12월 15일 - Korea Evangelist - Site Home - MSDN Blogs

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s