This past weekend, I needed to find a good fraction to approximate a given value. Along the way, I learned about continued fractions, and emerged with a neat, efficient algorithm. The only trouble is, I can’t figure out why it generates the correct answer!

The application

I have a data compression application in mind, based on linear interpolation between two points. Since Vena has already developed a library for unlimited-precision fractions, I figured I’d use that to pre-compute the slope between the points, and then I’d interpolate the “rise” by multiplying the “run” by the slope.

Our fraction class has two implementations: one called BigFraction using Java’s BigInteger for the numerator and denominator, and a much faster one called SmallFraction based on long for reasonably small fractions. For performance reasons, I’d like to use SmallFraction whenever possible.

Since the interpolation can tolerate small errors, there’s no point in using an enormously precise slope fraction. We can limit the numerator and denominator to, say, six digits; that way, even when we multiply by the “run” (which is limited to be a 32-bit int), the result still fits in a SmallFraction.

Now, how should we approximate the slope with a fraction having its numerator and denominator limited to six digits? Here’s one straightforward way:

• Take the numerator (or denominator, if that’s larger)
• If it’s under one million, stop
• Otherwise, scale it to 999,999, and scale the denominator (or numerator) accordingly

For example, with this technique, pi would become 999999/318309 which is off by about two parts per million.

Simple, right? But on the weekend, I started to wonder whether this was the best I could do.

Continued fractions

It turns out the problem I’m trying to solve is called the best rational approximation problem, which is closely related to a very cool concept called continued fractions. I don’t have the space in this blog post to describe continued fractions in detail, but the Wikipedia page goes a very good job of explaining them.

Here’s an example of a continued fraction representing pi:

This is a pretty cumbersome way to write a fraction, so instead, they write it using a special notation, with the whole number part followed by a semicolon, followed by a comma-separated list of the numbers you’d see in the denominator. For pi, we’d write:[3; 7, 15, 1, 292, 1, 1, 1, …]

The reason this interests me is that wherever you get tired of writing a continued fraction and decide to stop, you’ll have a new record for best approximation so far. For pi, the best integer approximation is 3. The next-best approximation^† is the continued fraction [3; 7], which is 3+1/7, or 22/7. The next-best is [3; 7, 15], which is 3+1/(7+1/15), or 333/106; followed by [3; 7, 15, 1] which is 355/113. (This last fraction has half as many digits as the 12-digit monster we computed before, yet it’s 20 times closer to the right value!)

In general, the continued fraction written [a; b, c, d, e] would correspond to this expression:a + 1/(b+1/(c+1/(d+1/(e))))

All those parentheses mean that, in effect, continued fractions are evaluated from right to left. This is going to pose a serious problem very shortly.

Since continued fractions are so good at approximating things, they lend themselves to a straightforward approximation algorithm, which goes something like this:

1. Initialize C to the continued fraction representation of the ideal slope.
2. Check whether the fraction representation of C has numerator and denominator under the limit. If so, return this fraction.
3. Else, chop off the last element of C and try again.

This is a pretty crummy algorithm. It has two major problems:

Problem 1: Step 1 starts by computing the whole continued-fraction representation of the given slope, even though we don’t need it all, which can be very inefficient; in fact, for irrational numbers, the continued fraction never terminates.

To solve this problem, we’d like our algorithm to interleave the generation of the continued fraction and the limit check. That leads us to a new algorithm that looks more like this:

1. Initialize F to the integer portion of the ideal slope, and C to the corresponding integer-only continued-fraction representation.
2. Compute the next term in C (as described on Wikipedia) and the corresponding fraction G.
3. If numerator or denominator of G is past the limit, we’ve gone too far. Return the previous best fraction F.
4. Otherwise, set F to G and repeat.

However, this doesn’t solve the second major problem with this approach.

Problem 2: Step 2 computes the fraction representation of a given continued fraction from scratch each time. Since this is an O(n) operation, the algorithm is O(n²).

This is a tough one. As long as we’re forced to evaluate continued fractions from right to left, it means that partial results from previous loop iterations don’t help us. If you want to compute 3+1/(7+1/(15+1/(1+1/(292)))), it doesn’t do you much good to know that the previous answer, 3+1/(7+1/(15+1/(1))), happens to be 355/113.

Or does it?

Let’s call it empirical mathematics

This is the point where I left the path of deductive reasoning, and started messing around with examples until I noticed a pattern.

Suppose we want to compute [0; a₁, a₂, …, a_n], and let’s call this C_n. The problem is the following: how do we compute C_n given a_n, and information about all previous C calculations?

Here’s the answer I found. Define two sequences N_n and D_n as follows:
N_n = a_nN_n-1 + N_n-2
D_n = a_nD_n-1 + D_n-2

These two equations look very much like a generalized form of the Fibonacci sequence, where F_n = F_n-1 + F_n-2. They’re simple and efficient to calculate.

But here’s the thing: Initialize these sequences properly, and lo and behold, you’ll find that C_n = N_n/D_n. I subsequently found this formula on Wolfram MathWorld, so at least I’m able to confirm that it always works.

What we’ve got here is a simple, efficient way to compute each C_n in constant time and constant space^‡. The algorithm now looks like this:

1. Initialize N₀ = 0, D₀ = 1, N₁ = 1, D₁ = a₁
2. Compute the next continued fraction term, a_n, as described on Wikipedia (starting with n = 2)
3. Set N₂ = a_nN₁ + N₀ and D₂ = a_nD₁ + D₀
4. If N₂ or D₂ is past the limit, we’ve gone too far; return N₁ / D₁
5. Otherwise, shift everything down (so D₂ becomes D₁ and so forth), and repeat

These calculations do a fixed amount of work each time around the loop; they operate entirely on small integers that can fit in machine registers; and they’re simple and fast. It’s a great algorithm, and I’m altogether quite impressed with myself for having (re-)invented it.

Now, could someone please tell me how it works?

†There are actually a few other pretty good approximations in between these “best” ones, but we’ll ignore those for simplicity. It’s all on the Wikipedia page.

‡Assuming a machine model in which integers occupy a fixed amount of space, and multiplies and divides occur in constant time. This is a pretty loose interpretation of asymptotic complexity, but we can get away with it here if we assume that the algorithm to terminates while all the numbers are small enough to make these assumptions correct (which it provably always does).

• Category: Vena Engineering
• Tag: algorithms, math