ALGORITHM ALLEY

Continued Fractions Versus Farey Fractions

Andreas Bender

Andreas is a Ph.D. student in the University of Cambridge's Department of Pure Mathematics and Mathematical Statistics. He can be contacted at A.O.Bender@ pmms.cam.ac.uk.

Introduction

Speed and accuracy are common concerns in the study of algorithms. In time-critical code, it can even be difficult to do accurate arithmetic. Say, for example, you need to multiply by in a fixed-point computation. One obvious approach is to convert into a fixed-point number, and do a single multiplication. Unfortunately, accuracy requires working with large intermediate values. It's often better to multiply by a small value (to avoid overflow) and divide by another small value. If you choose your fraction carefully, you can actually obtain better precision than with the fixed-point approach.

In his October 1995 "Algorithm Alley," Louis Plebani examined one way of finding such fractions using Farey series. This month, Andreas Bender takes a closer look at the Farey series and presents an alternative approach based on continued fractions.

--Tim Kientzle

Writing a number such as ".12345" as a fraction is entirely straightforward--=12345/100000. However, what fraction best approximates if its numerator or denominator may not exceed a given size H? How small can the numerator or denominator of the approximating fraction be if we can tolerate an error of at most .

Careful readers of this column will remember that Louis Plebani addressed these questions in his article "Common-Fraction Approximation of Real Numbers" (DDJ, October 1995). Instead of the Farey fractions he used, however, I'll approach the problem using continued fractions, then compare the merits of the two methods. While I hope that this article is reasonably self contained, it will certainly be easier to follow if you have read Plebani's article as well.

What is a continued fraction? Figure 1 provides a concise definition. To build a continued-fraction approximation to a rational number, we apply the Euclidean algorithm to the numerator and denominator of =r/s; see Figure 2. Since <1, q₀ is equal to 0. The numbers q_i are called the "partial quotients." Truncating the continued fraction of Figure 1 at q_i leaves the "complete quotient" n_i/d_i. For example, the complete quotients of ¹²/₁₇ are ⁰/₁, ¹/₁, ²/₃, ⁵/₇, and ¹²/₁₇. The theory of Diophantine approximation produces the results:

Every complete quotient n_i/d_i is the best approximation to f with denominator d_i.
The error |(n_i/d_i)-| is less than, or equal to, 1/d_i².

This suggests the following procedure: Calculate one complete quotient after the other and test every one of them for our condition on approximation accuracy or size of denominator until it is satisfied.

The Farey approximation, as described by Plebani, is calculated by subdividing an interval. The interval from a/b to a'/b' is subdivided at (a+a')/(b+b'); for the start interval [0,1] this would, of course, be 1/2. You then check in which of the two newly generated subintervals lies. A median can now be inserted into this interval and this process iterated until some endpoint satisfies the desired condition.

The important facts about Farey series are as follows:

The n^th Farey series contains all fractions between 0 and 1 with denominator at most n.
There is a u/v in this series for which |u/v - |1/((n+1)v) holds.
This u/v, calculated by the aforementioned procedure, is the best approximation of by a fraction with denominator at most n.

These two approaches have an important connection: Each subsequent complete quotient in the continued-fraction calculation gives the same result as computing q_i medians (see R.L. Graham et al, Concrete Mathematics, Addison-Wesley, 1989). As a result, the continued-fraction approach gives more accurate results more quickly, but the Farey series approach is more sensitive to bounds on the size of the denominator.

It's possible to mix these two approaches, starting with a continued-fraction estimate, then extending the result to the optimal Farey approximation. To do this, I need to find the successor (or predecessor) of n_l/d_lin the Farey series of order d_l.Start with the complete quotient n_l/d_land assume that > n_l/d_l (for the other case, only the sign gcd(r,s) needs to be reversed). By tracing the steps of the Euclidean algorithm, you can solve the equation d_lx-n_ly=gcd(r,s). Now divide the whole equation by the gcd(r,s) and, henceforth, assume the gcd to be one. If (x₀,y₀) is a solution, then (x₀+jn_l,y₀+jd_l) also is one for every integer j, so we can pick a j such that 0<y₀+jd_l n_l and this j gives us the successor of n_l/d_l in the Farey series of order d_l. Note that the most computation-intensive part of this procedure--finding (x₀,y₀)--amounts to just finishing the Euclidean algorithm for r and s and tracking the partial quotients. That is the algorithm used for computing the complete quotients as well.

The larger the q_i, the more interesting the question of grouping these q_i steps together becomes. Given an arbitrarily chosen , what are the likely sizes of the q_i? A thorough discussion of this question entails delicate arithmetic estimates and an examination of the underlying assumptions (see A.M. Rockett et al., "Continued Fractions," World Scientific, 1992).

A crude summary is that q_i is "not bounded, but not too big too often either." In The Art of Computer Programming, Vol.2: Seminumerical Algorithms (Addison-|Wesley, 1969), Donald Knuth notes two interesting facts:

q_i is equal to a with probability P(a)=log₂ ((a+1)²/(a+1)²-1). The first few values are P(1)=0.415, P(2)=0.17, P(3)=0.093, P(4)=0.059.
If the denominator is equal to s, we can expect the Euclidean algorithm to take about 1.94 log₁₀s steps, where I disregard certain correction terms.