blog

The modern theory of continued fractions comes from Christiaan Huygens, a Dutch physicist who invented the pendulum clock. Continued fractions turn out to be an especially elegant way of finding rational approximations of a number; this enabled him to design clocks with small gears that nonetheless provided the desired degree of accuracy.

Several hundred years later, continued fractions are still the best way to compute rational approximants. But as we no longer compute them by hand, there is likewise no reason to compute them with while loops; continued fractions lend themselves to an elegant treatment using apomorphisms.

Every real number has an expansion of the form shown below. Hence, a continued fraction is completely determined by a list of positive integers; we can use [Integer] to represent it in Haskell.

$$[a_o;a_1,a_2,\ldots] = a_o + \frac{1}{a_1 + \frac{1}{a_2 + \cdots}}$$

We start with the prolegomena, and then the continuedFraction function we will use for conversions.

import           Data.Functor.Foldable
isInteger :: (RealFrac a) => a -> Bool
isInteger = idem (realToFrac . floor)
    where idem = ((==) <*>)
continuedFraction :: (RealFrac a, Integral b) => a -> [b]
continuedFraction = apo coalgebra
    where coalgebra x
              | isInteger x = go $ Left []
              | otherwise   = go $ Right alpha
                  where alpha = 1 / (x - realToFrac (floor x))
                        go    = Cons (floor x)

This is mostly straightforward, despite the abuse of pointfree notation in idem. continuedFraction could have been written via pattern matching as well, but it is more elegant and certainly more instructive this way. Though apormophisms are perhaps more well-known than Elgot algebras, concrete examples can nonetheless be valuable in this area.

Next, we need to write a function to convert continued fractions back to a form we understand. We'll just use vanilla pattern matching this time, but it works fairly nicely here.

collapseFraction :: (Integral a) => [Integer] -> Ratio a
collapseFraction [x]    = fromIntegral x % 1
collapseFraction (x:xs) = (fromIntegral x % 1) + 1 / collapseFraction xs

We now make use of the (mathematical) fact that the denominators of the convergents of a continued fraction form a strictly increasing sequence, and the fact that subsequent convergents always have a smaller error. Hence we can simply consider all rational approximations, stopping when the denominator becomes too large, and this will produce the best approximation, as desired.

The rest falls out naturally. Note that continuedFraction will in general return an infinite list, so any time we use it, we have to compose it with take.

convergent :: (RealFrac a, Integral b) => a -> Integer -> Ratio b
convergent x n = collapseFraction $ take n (continuedFraction x)
approximate :: (RealFrac a, Integral b) => a -> b -> Ratio b
approximate x d = last . takeWhile ((<= d) . denominator) $ fmap (convergent x) [1..]

At this point, we have exactly what we want. A simple approximate pi 100 will yield the best rational approximation of π with a denominator less than 100, a familiar 22 % 7.