Using the Fractional Type that we introduced in a previous blogpost, we build a module for continued fractions in Haskell in this entry. As usual, the code is licensed BSD and is authored by myself. Nonetheless, I would like to thank Cale Gibbard for his invaluable input for the debugging of certain parts of it. Since this is quite a big blog post, I’ve decided to post it already partially finished, as I have a short attention span but wanted to get the code out there. So be warned that this is WIP, and is not finished yet. I will remove this notice when it is.
Originally, I read about continued fractions on “Arithmetic with Continued Fractions” by Dominus. Intrigued by them, I digged a bit further, but sadly was unable to find much more than the original paper translated in plain html with rather old and hard to read notation. What intrigued me about this notation is that one could get real numbers to an unlimited precision in a lazy language like Haskell. After all, using laziness, we only compute what we really need. So once we have some continued fraction, it is easy to imagine only getting N parts of it to get a value within some eps of the final value. This is where this module originally originated. Another interesting aspect of continued fractions is that they can be used to represent both rational and irrational numbers to any precision. Additionally, a nice feature of them is that all rational numbers have a finite representation with them.
After having completely implemented the continued fractions, I noticed that it was actually really hard to deal with irrational numbers. The current implementation of continued fractions simply does not support them. That being said, I think it is a simple extension of this implementation that can support them rather well. I will focus on this and try to write it down here after detailing the implementation as it is.
Continued Fraction Representations
There are various representations for continued fractions, as can be seen in the list below. All of them allow for representing any number, however some of them are more convenient than others.
- a0 + b0a1 + b1a2 + b2…
- a0 + 1a1 + 1a2 + 1…
The second representation is often named regular continued fractions. The advantage of the prior representation is that it’s easier to represent certain irrational numbers such as pi:
- π4 = 11+13+45+97+16/…
- pi4 : a0 = 0, b0 = 1, ∀ i > 0.(ai, bi) = (2·i−1,i2)
That being said, the advantage of the latter representation is that all numbers are uniquely representable (well not quite exactly, but we’ll look at that later).
Continued Fractions Arithmetic
TODO
Continued Fractions Implementation
Using the custom definition of Fraction from one of my previous blogposts we can now go ahead and define the actual code for continued fractions. First we start with the header of the file, with the license, some specific GHC flags and some imports:
> {-# LANGUAGE BangPatterns #-}
> ------------------------------------------------------------------------------
> -- CF.hs - Continued Fractions
> --
> -- Authors: Christophe Poucet - <christophe (dot) poucet (at) gmail (dot) com>
> -- License: BSD
> -- Created: 01/22/2008 06:59:59 AM CET, last modified 01/25/2008
> --
> -- Copyright 2008 © Christophe Poucet. All Rights Reserved.
> -- Thanks to Cale Gibbard for his invaluable input
> ------------------------------------------------------------------------------
> module CF where -- (CF, approximate, approximations) where
> import Fraction
> import Data.List
> import Data.Ratio
The first important data type, obviously, is the continued fraction, which is a list of digits. Additionally, we want two structures to store the transformation ‘matrices’ for the arithmetic of one or two continued fractions. These are T4 and T8 respectively.
> --
> -- * Data definitions
> --> newtype CF a = CF [a]
> deriving (Eq, Show)
>
> data T4 a = T4 !a !a !a !a
> deriving (Eq)
>
> data T8 a = T8 !a !a !a !a !a !a !a !a
> deriving (Eq)
Simple function to turn a continued fraction into an exact value. The code from this comes from “Functional Pearl – Enumerating Rationals” by Jeremy Gibbons et al. It is much more efficient than the original solution I had.
> --
> -- | exact : Calculates the exact value of a continued fraction
> --
> exact :: (Integral a, Fractional b) => CF a -> b
> exact (CF ps) = a/b
> where (a, b) = foldr op (1,0) . map fromIntegral $ ps
> op m (n,d) = (m*n+d, n)
Checks whether a continued fraction is infinity. Notice that the for the case that any digits are negative, the logic is cruder. Then again, we do not allow negative digits, and this is simply a remnant from some deep debugging.
> --
> -- | isInfinity : Checks whether a CF is infinity
> --
> isInfinity (CF xs) = check xs || (any (<0) $ drop 1 xs)
> where check [] = True
> check (_:0:xs) = check xs
> check _ = False
TODO
> ------------------------------------------------------------------------------
> -- * Core Functions
> ------------------------------------------------------------------------------
> ------------------------------------------------------------------------------
> --
> -- | morph : allows for scaling and addition of normal numbers to a CF
> --
> -- a + bx
> -- (T4 a b c d) : z = --------
> -- c + dx
> --
> ------------------------------------------------------------------------------
> morph :: (Integral a) => T4 a -> CF a -> CF a
> morph (T4 a b c d) (CF xs) = CF $ worker a b c d xs
> where worker !a !b !c !d xs =
> --------------------------------------------------------------------
> case () of > _ | infinity -> []
> ------------------------------------------------------------------
> _ | fixed -> z:worker a' b' c' d' xs
> where a' = c
> b' = d> c' = (a-c*z)
> d' = (b-d*z)
> ------------------------------------------------------------------
> _ | not ended -> worker a' b' c' d' (tail xs)
> where a' = b> b' = (a+b*x)
> c' = d> d' = (c+d*x)
> ------------------------------------------------------------------
> _ | otherwise -> worker a' b' c' d' xs
> where a' = b
> b' = b
> c' = d
> d' = d
> --------------------------------------------------------------------
> where qs@[q1, q2] = case () of
> _ | ended -> [a :~: c, b :~: d]
> _ | otherwise -> [(a + b*x ) :~: (c + d*x) ,
> (a + b*(x+1)) :~: (c + d*(x+1)) ]
> infinity = c == 0 && d == 0
> ended = null xs
> all_equal = equal . filter (not . isNan) $ qs
> -- Note, lossy-fraction's equality is -NOT- transitive
> fixed = not infinity && not vanishing && all_equal
> -- Vanishing: We do not cross 0 in our quadrant, subject to 'equal qs'
> vanishing = not ended && signum c /= (signum d * signum x)
> x = head xs
> z = toIntegral . head . filter (not . isNan) $ qs
TODO
> ------------------------------------------------------------------------------
> --
> -- | arithmetic : allows for arithmetic operators between two CFs'
> --
> -- a + bx + cy + dxy
> -- (T8 a b c d e f g h) : z = -------------------
> -- e + fx + gy + hxy
> --
> ------------------------------------------------------------------------------
> arithmetic :: (Integral a) => T8 a -> CF a -> CF a -> CF a
> arithmetic (T8 a b c d e f g h) (CF xs) (CF ys) =
> CF $ worker a b c d e f g h xs ys
> where worker :: (Integral a) =>
> a -> a -> a -> a ->
> a -> a -> a -> a ->
> [a] -> [a] -> [a]
> worker !a !b !c !d !e !f !g !h xs ys =
> --------------------------------------------------------------------
> case () of
> _ | infinity -> []
> ------------------------------------------------------------------
> -- We have an integer value that is the same for the entire quadrant
> _ | fixed ->
> z:worker a' b' c' d' e' f' g' h' xs ys
> where {
> a' = e ; b' = f ;
> c' = g ; d' = h ;
> e' = (a-e*z); f' = (b-f*z);
> g' = (c-g*z); h' = (d-h*z)
> }
> ------------------------------------------------------------------
> -- x /= inf || y /= inf
> _ | not ended ->
> --------------------------------------------------------------
> case () of
> ------------------------------------------------------------
> -- X shows largest variance
> _ | choose_x ->
> case () of
> ------------------------------------------------------
> -- x == inf, change variance along x to 0
> _ | ended_x ->
> worker a' b' c' d' e' f' g' h' xs ys
> where {
> a' = 0 ; b' = b ;
> c' = 0 ; d' = d ;
> e' = 0 ; f' = f ;
> g' = 0 ; h' = h
> }
> ------------------------------------------------------
> -- x /= inf
> _ | otherwise ->
> worker a' b' c' d' e' f' g' h' (tail xs) ys
> where {
> a' = b ; b' = (a+b*x);
> c' = d ; d' = (c+d*x);
> e' = f ; f' = (e+f*x);
> g' = h ; h' = (g+h*x)
> }
> ------------------------------------------------------------
> -- Y shows largest variance
> _ | otherwise ->
> case () of
> ------------------------------------------------------
> -- y == inf, change variance along y to 0
> _ | ended_y ->
> worker a' b' c' d' e' f' g' h' xs ys
> where {
> a' = 0 ; b' = 0 ;
> c' = c ; d' = d ;
> e' = 0 ; f' = 0 ;
> g' = g ; h' = h
> }
> ------------------------------------------------------
> -- y /= inf
> _ | otherwise ->
> worker a' b' c' d' e' f' g' h' xs (tail ys)
> where {
> a' = c ; b' = d ;
> c' = (a+c*y); d' = (b+d*y);
> e' = g ; f' = h ;
> g' = (e+g*y); h' = (f+h*y)
> }
> --------------------------------------------------------------
> where choose_x | numNans >= 3 = not ended_x
> | isNan d_x && isNan d_y =
> d_xy_x > d_xy_y
> | isNan d_x && isNan d_xy_y =
> d_xy_x > d_y
> | isNan d_xy_x && isNan d_y =
> d_x > d_xy_y
> | isNan d_xy_x && isNan d_xy_y =
> d_x > d_y
> | numNans >= 2 = not ended_x
> | isNan d_x =
> d_xy_x > max d_y d_xy_y
> | isNan d_xy_x =
> d_x > max d_y d_xy_y
> | isNan d_y =
> max d_x d_xy_x > d_xy_y
> | isNan d_xy_y =
> max d_x d_xy_x > d_y
> | otherwise =
> max d_x d_xy_x > max d_y d_xy_y
> ds = [d_x, d_y, d_xy_x, d_xy_y]
> numNans = length . filter isNan $ ds
> d_x = abs $ q2 - q1
> d_y = abs $ q3 - q1
> d_xy_x = abs $ q4 - q3
> d_xy_y = abs $ q4 - q2
> ------------------------------------------------------------------
> -- x == inf && y == inf
> _ | otherwise ->
> worker a' b' c' d' e' f' g' h' xs ys
> where {
> a' = d ; b' = d ;
> c' = d ; d' = d ;
> e' = h ; f' = h ;
> g' = h ; h' = h
> }
> --------------------------------------------------------------------
> where qs@[q1, q2, q3, q4] = case () of
> _ | ended -> [a :~: e, b :~: f, c :~: g, d :~: h]
> _ | ended_x -> [(a + c*y ) :~: (e + g*y ) ,
> (b + d*y ) :~: (f + h*y ) ,
> (a + c*(y+1)) :~: (e + g*(y+1)) ,
> (b + d*(y+1)) :~: (f + h*(y+1)) ]
> _ | ended_y -> [(a + b*x ) :~: (e + f*x ) ,
> (a + b*(x+1)) :~: (e + f*(x+1)) ,
> (c + d*x ) :~: (g + h*x ) ,
> (c + d*(x+1)) :~: (g + h*(x+1)) ]
> _ | otherwise -> [(a + b*x + c*y + d*x*y )
> :~: (e + f*x + g*y + h*x*y ) ,
> (a + b*(x+1) + c*y + d*(x+1)*y )
> :~: (e + f*(x+1) + g*y + h*(x+1)*y ) ,
> (a + b*x + c*(y+1) + d*x*(y+1) )
> :~: (e + f*x + g*(y+1) + h*x*(y+1) ) ,
> (a + b*(x+1) + c*(y+1) + d*(x+1)*(y+1))
> :~: (e + f*(x+1) + g*(y+1) + h*(x+1)*(y+1)) ]
>
> infinity = e == 0 && f == 0 && g == 0 && h == 0
> all_equal = equal . filter (not . isNan) $ qs
> -- Note, lossy-fraction's equality is -NOT- transitive
> ended_x = null xs
> ended_y = null ys
> ended = ended_x && ended_y
> -- Note, lossy-fraction's equality is -NOT- transitive
> fixed = not infinity && not vanishing && all_equal
> -- Vanishing: We do not cross 0 in our quadrant, subject to 'equal qs'
> vanishing = case () of
> --------------------------------------------------------------
> _ | not ended_x && not ended_y ->
> -- -e - f*x
> -- y(x) = ----------
> -- g + h*x
> --
> within y (-e - f*x) (g + h*x) ||
> withinS y (-e - f*(x+1)) (g + h*(x+1)) ||
> -- -e - g*y
> -- x(y) = ----------
> -- f + h*y
> --
> within x (-e - g*y) (f + h*y) ||
> withinS x (-e - g*(y+1)) (f + h*(y+1))
> --------------------------------------------------------------
> _ | not ended_x ->
> signum g /= (signum h * signum x)
> --------------------------------------------------------------
> _ | not ended_y ->
> signum f /= (signum h * signum y)
> --------------------------------------------------------------
> _ | otherwise ->
> False
> x = head xs
> y = head ys
> z = toIntegral . head . filter (not . isNan) $ qs
> --
> -- ** Comparison
> --
> instance (Ord a) => Ord (CF a) where
> compare (CF ps) (CF qs) = comp ps qs
> where comp [] [] = EQ
> comp (x:xs) [] = LT
> comp [] (y:ys) = GT
> comp (x:xs) (y:ys) = if c == EQ then comp ys xs else c
> where c = compare x y
> --
> -- ** Standard Arithmetic operators
> --
> instance (Integral a) => Num (CF a) where
> x + y = arithmetic (T8 0 1 1 0 1 0 0 0) x y
> x * y = arithmetic (T8 0 0 0 1 1 0 0 0) x y
> x - y = arithmetic (T8 0 1 (-1) 0 1 0 0 0) x y
> negate (CF ps) = CF $ worker ps
> where worker [] = []
> worker [a] = [-a]
> worker [a,2] = [-a-1,2]
> worker (a:1:b:r) = (-a-1):(b+1):r
> worker (a:b:r) = (-a-1):1:(b-1):r
> abs (CF []) = CF []
> abs cx@(CF (x:xs)) | x < 0 = negate cx
> | otherwise = cx
> signum (CF []) = fromInteger 1
> signum (CF (x:xs)) | x < 0 = fromInteger (-1)
> | x == 0 &&
> null xs = fromInteger 0
> | otherwise = fromInteger 1
> fromInteger x = CF [fromInteger x]
> --
> -- ** Enumeration
> --
> instance (Integral a) => Enum (CF a) where
> succ x = morph (T4 1 1 1 0) x
> pred x = morph (T4 (-1) 1 1 0) x
> toEnum x = fromInteger . toInteger $ x
> fromEnum (CF []) = error "fromEnum: The CF is infinity"
> fromEnum (CF (x:xs)) = fromIntegral x
> enumFrom x = iterate succ x
> enumFromThen x y = iterate ((y-x)+) x
> enumFromTo x y = takeWhile (<= y) $ enumFrom x
> enumFromThenTo x y z = takeWhile p (enumFromThen x y)
> where p | y >= x = (<= z)
> | otherwise = (>= z)
> --
> -- ** Fractional
> --
> instance (Integral a) => Fractional (CF a) where
> x / y = arithmetic (T8 0 1 0 0 0 0 1 0) x y
> recip x = morph (T4 1 0 0 1) x
> fromRational x = CF $ loop (numerator x) (denominator x)
> where loop n b | b == 0 = []
> | otherwise = let (d,m) = n `divMod` b
> in fromIntegral d:loop b m
> instance (Integral a) => Real (CF a) where
> toRational cf = exact cf
> ------------------------------------------------------------------------------
> -- * Helper functions
> ------------------------------------------------------------------------------
> none :: (a -> Bool) -> [a] -> Bool
> none f l = not . any f $ l
> {-# INLINE none #-}
>
> equal :: (Eq a) => [a] -> Bool
> equal (x:xs) = all (x ==) xs
> equal _ = True> {-# INLINE equal #-}
>
> -- within z a b == z <= a/b < z+1
> -- Could be rephrased as a CF function: a :~: b == z :~: 1
> within z a b | b == 0 = a == 0
> | b > 0 = b*z <= a && a < b*(z+1)
> | otherwise = b*z >= a && a > b*(z+1)
> {-# INLINE within #-}
>
> -- withinS z a b == z < a/b < z+1
> -- Could be rephrased as a CF function: a :~: b == z :~: 1 && a /= z*b
> withinS z a b | b == 0 = a == 0
> | b > 0 = b*z < a && a < b*(z+1)
> | otherwise = b*z > a && a > b*(z+1)
> {-# INLINE withinS #-}
A simple test to see whether a continued fraction is normalized (something very crucial during the debugging of two CF arithmetic) is simply to check what happens when you multiply it by 1.
> --
> -- | isNormalized : Whether a CF is in the proper normal form
> --
> isNormalized cf@(CF xs) = xs == ys
> where (CF ys) = morph (T4 0 1 1 0) cf
Finally, we simply define some show instances for the T4 and T8 data structure, to ease algorithmic debugging (which was a great must during the development of this library).
> instance Show a => Show (T4 a) where
> show (T4 a b c d) = unlines [unwords ["|/", as, bs, "\\"],
> unwords ["|\\", cs, ds, "/"]]
> where s = map show [a,b,c,d]
> len = maximum . map length $ s
> [as,bs,cs,ds] = map (\x -> take (len - length x) spaces ++ x) s
> spaces = cycle " "
>
> instance Show a => Show (T8 a) where
> show (T8 a b c d e f g h) = unlines [unwords ["|/", as, bs, cs, ds, "\\"],
> unwords ["|\\", es, fs, gs, hs, "/"]]
> where s = map show [a,b,c,d,e,f,g,h]
> len = maximum . map length $ s
> [as,bs,cs,ds,es,fs,gs,hs] = map (\x -> take (len - length x) spaces ++ x) s
> spaces = cycle " "
Some final thoughts
Now as mentioned, I’d explain some issues regarding the use of this for irrationals and a possible solution. After implementing this system, I was obviously curious to try this out on an irrational. So I tried multiplying the CF for √2 with itself. And after seeing it not converge at all, I kept wondering why. Now I know, and it also shows a way to fix it so that it will work. Due to the representation I use for continued fractions, I am guaranteed that the integer part will always be stable. The problem is that I can not have negative additions, only positive ones. This means that if I were able to get the stable part ‘2′ early, then the remaining infinite digits of the CF for √2 would have to add up to 0. One possibility for ‘fixing’ this is by changing the measure of stability. If we allow negative digits, we could choose to stabilize a digit as long as the different quotients are all within some ε from a given integer, and not require them all to be falling on the same side. This would, of course, mean, that it is possible to have negative digits further in the stream. Whether we want this or not is a big question, since they open up a whole other can of worms regarding stability of the number and such, perhaps a more mathematically-inclined person could build further upon this idea.
