cf 0.1 → 0.2

raw patch · 5 files changed

+725/−180 lines, 5 filesdep +QuickCheckdep +cfdep +test-frameworkdep ~basePVP ok

version bump matches the API change (PVP)

Dependencies added: QuickCheck, cf, test-framework, test-framework-quickcheck2, test-framework-th

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• cf.cabal +20/−4
• src/Math/ContinuedFraction.hs +293/−176
• src/Math/ContinuedFraction/Interval.hs +163/−0
• src/Math/ContinuedFraction/Simple.hs +201/−0
• tests/Tests.hs +48/−0

@@ -1,6 +1,6 @@ name:                cf-version:             0.1-synopsis:            Infinite precision arithmetic using continued fractions+version:             0.2+synopsis:            Exact real arithmetic using continued fractions license:             MIT    license-file:        LICENSE author:              Mitchell Riley@@ -15,7 +15,23 @@   location: git://github.com/mvr/cf.git  library-  hs-source-dirs:      src-  exposed-modules:     Math.ContinuedFraction+  hs-source-dirs:  src+  exposed-modules: Math.ContinuedFraction,+                   Math.ContinuedFraction.Simple,+                   Math.ContinuedFraction.Interval   build-depends:       base >=4.7 && <4.8   default-language:    Haskell2010++test-suite tests+  type: exitcode-stdio-1.0+  main-is: Tests.hs+  default-language:    Haskell2010+  hs-source-dirs:+    tests+  build-depends:+    base,+    cf,+    QuickCheck                 >= 2.4,+    test-framework             >= 0.6,+    test-framework-quickcheck2 >= 0.2,+    test-framework-th          >= 0.2

@@ -1,216 +1,333 @@-module Math.ContinuedFraction-  (-    CF,-    convergents,-    digits,-    showCF,-    sqrt2,-    exp1-  ) where+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE TypeFamilies #-}+module Math.ContinuedFraction (+  CF,+  cfString,+  cfcf+) where +import Data.Maybe (catMaybes, mapMaybe) import Data.Ratio -newtype CF = CF [Integer]+import Math.ContinuedFraction.Interval --- | Produce a list of rational approximations to a number-convergents :: CF -> [Rational]-convergents (CF cf) = go 0 1 1 0 cf-  where go p q p' q' (a:as) = (newp % newq) : go p' q' newp newq as-          where newp = a * p' + p-                newq = a * q' + q-        go _ _ _ _ [] = []+newtype CF' a = CF [a]+type CF = CF' Integer --- | Produce a list of digits in a given base-digits :: Integer -> CF -> [Integer]-digits base (CF cf) = go 0 1 1 0 cf-  where go 0 _ 0 _ _        = []-        go _ _ p' q' []     = go p' q' p' q' [0]-        go p q p' q' (a:as) = case digit p q p' q' of-                                Just d -> d : go (base * (p - d * q)) q (base * (p' - d * q')) q' (a:as)-                                Nothing -> go p' q' (a * p' + p) (a * q' + q) as-        digit p q p' q' = if q' /= 0 && q /= 0 && p `quot` q == p' `quot` q' then-                            Just $ p `quot` q-                          else-                            Nothing+type Hom a = (a, a,+              a, a) --- | Produce a decimal representation of a number-showCF :: CF -> String-showCF cf | cf < 0 = "-" ++ show (-cf)-showCF (CF [i])   = show i-showCF (CF (i:r)) = show i ++ "." ++ decimalDigits-  where decimalDigits = concatMap show $ tail $ digits 10 (CF (0:r))+type Bihom a = (a, a, a, a,+                a, a, a, a) --- Should make this cleverer-instance Show CF where-  show = take 15 . showCF+class (Fractional (FractionField a)) => HasFractionField a where+  type FractionField a :: *+  insert :: a -> FractionField a+  extract :: FractionField a -> (a, a) -safeHead :: [a] -> Maybe a-safeHead (x:_) = Just x-safeHead [] = Nothing+instance HasFractionField Integer where+  type FractionField Integer = Rational+  insert = fromInteger+  extract r = (numerator r, denominator r) -safeRest :: [a] -> [a]-safeRest (_:xs) = xs-safeRest [] = []+instance HasFractionField Rational where+  type FractionField Rational = Rational+  insert = id+  extract r = (numerator r % 1, denominator r % 1) --- The coefficients of the homographic function (a + bx) / (c+dx)-type Hom = (Integer, Integer,-            Integer, Integer)+instance HasFractionField CF where+  type FractionField CF = CF+  insert = id+  extract r = (r, 1) --- Possibly output a term and return the simplified hom-emit :: Hom -> Maybe (Hom, Integer)-emit (a, b,-      c, d) = if c /= 0 && d /= 0 && r == s then-                Just ((c,       d,-                       a - c*r, b-d*r), r)-              else-                Nothing-  where r = a `quot` c-        s = b `quot` d+homEmit :: Num a => Hom a -> a -> Hom a+homEmit (n0, n1,+         d0, d1) x = (d0,        d1,+                      n0 - d0*x, n1 - d1*x) --- Absorb the next term-ingest :: Hom -> Maybe Integer -> Hom-ingest (a, b,-        c, d) (Just p) = (b, a+b*p,-                          d, c+d*p)-ingest (_a, b,-        _c, d) Nothing  = (b, b,-                           d, d)+homAbsorb :: Num a => Hom a -> a -> Hom a+homAbsorb (n0, n1,+           d0, d1) x = (n0*x + n1, n0,+                        d0*x + d1, d0) --- Apply a hom to a continued fraction-hom' :: Hom -> [Integer] -> [Integer]-hom' (0, 0,-      _, _) _ = [0]-hom' (_, _,-      0, 0) _ = []-hom' h x = case emit h of-           Just (next, d) -> d : hom' next x-           Nothing -> hom' (ingest h (safeHead x)) (safeRest x)+det :: Num a => Hom a -> a+det (n0, n1,+     d0, d1) = n0 * d1 - n1 * d0 -hom :: Hom -> CF -> CF-hom h (CF x) = CF $ hom' h x+homEval :: (Num a, HasFractionField a, Eq (FractionField a)) => Hom a -> Extended (FractionField a) -> Extended (FractionField a)+homEval (n0, n1,+         d0, d1) (Finite q) | denom /= 0 = Finite $ num / denom+                            | num == 0 = error "0/0 in homQ"+                            | otherwise = Infinity+  where num   = insert n0 * q + insert n1+        denom = insert d0 * q + insert d1+homEval (n0, _n1,+         d0, _d1) Infinity = Finite $ insert n0 / insert d0 --- The coefficients of the bihomographic function (a + bx + cy + dxy) / (e + fx + gy + hxy)-type Bihom = (Integer, Integer, Integer, Integer,-              Integer, Integer, Integer, Integer)+constantFor :: (Eq a, Num a, HasFractionField a) => Hom a -> Extended (FractionField a)+constantFor (_, _,+             0, 0) = Infinity+constantFor (0, 0,+             0, _) = Finite 0+constantFor (0, 0,+             _, 0) = Finite 0+constantFor (a, 0,+             b, 0) = Finite (insert a / insert b)+constantFor (_, a,+             _, b) = Finite (insert a / insert b) --- Possibly output a term and return the simplified bihom-biemit :: Bihom -> Maybe (Bihom, Integer)-biemit (a, b, c, d,-        e, f, g, h) = if e /= 0 && f /= 0 && g /= 0 && h /= 0 && ratiosAgree then-                      Just ((e,     f,     g,     h ,-                             a-e*r, b-f*r, c-g*r, d-h*r), r)+boundHom :: (Ord a, Num a, HasFractionField a, Eq (FractionField a)) => Hom a -> Interval (FractionField a) -> Interval (FractionField a)+boundHom h (Interval i s) | det h > 0 = Interval i' s'+                          | det h < 0 = Interval s' i'+                          | otherwise = Interval c c+  where i' = homEval h i+        s' = homEval h s+        c = constantFor h++primitiveBound :: forall a. (Ord a, Num a, HasFractionField a) => a -> Interval (FractionField a)+primitiveBound n | abs n < 1 = Interval (Finite $ insert bot) (Finite $ insert top)+  where bot = (-2) :: a+        top = 2 :: a+primitiveBound n = Interval (Finite $ an - 0.5) (Finite $ 0.5 - an)+  where an = insert $ abs n++-- TODO: just take the rational answer from the hom+nthPrimitiveBounds :: (Ord a, Num a, HasFractionField a, Eq (FractionField a)) =>+                       CF' a -> [Interval (FractionField a)]+nthPrimitiveBounds (CF cf) = zipWith boundHom homs (map primitiveBound cf) ++ repeat (Interval ev ev)+  where homs = scanl homAbsorb (1,0,0,1) cf+        ev = evaluate (CF cf)++evaluate :: (HasFractionField a, Eq (FractionField a)) => CF' a -> Extended (FractionField a)+evaluate (CF []) = Infinity+evaluate (CF [c]) = Finite $ insert c+evaluate (CF (c:cs)) = case next of+                        (Finite 0) -> Infinity+                        Infinity   -> Finite $ insert c+                        (Finite r) -> Finite $ insert c + recip r+  where next = evaluate (CF cs)++valueToCF :: RealFrac a => a -> CF+valueToCF r = if rest == 0 then+                CF [d]+              else+                let (CF ds)  = valueToCF (recip rest) in CF (d:ds)+  where (d, rest) = properFraction r++intervalThin :: (RealFrac a) => Interval a -> Bool+intervalThin (Interval Infinity    Infinity)  = False+intervalThin (Interval Infinity   (Finite _)) = False+intervalThin (Interval (Finite _)  Infinity)  = False+intervalThin (Interval (Finite i) (Finite s)) = abs z > 3 || abs (zi - zs) < 2+  where zi = round i+        zs = round s+        z  = if abs zs < abs zi then zs else zi++euclideanPart :: (RealFrac a, Integral b) => Interval a -> Maybe b+euclideanPart (Interval Infinity    Infinity)  = undefined+euclideanPart (Interval Infinity   (Finite b)) = Just $ floor b+euclideanPart (Interval (Finite a)  Infinity)  = Just $ ceiling a+euclideanPart i@(Interval (Finite a) (Finite b))+  | 0 `elementOf` i && not subsetZero = Nothing+  | zi /= 0 && zs /= 0 = Just z+  | subsetZero = Just 0+  | otherwise = Nothing+    where zi = round a+          zs = round b+          z  = if abs zs < abs zi then zs else zi+          subsetZero = i `subset` Interval (Finite (-2)) (Finite 2)++existsEmittable :: RealFrac a => Interval a -> Maybe Integer+existsEmittable i = if intervalThin i then+                      euclideanPart i                     else                       Nothing-  where r = a `quot` e-        ratiosAgree = r == b `quot` f && r == c `quot` g && r == d `quot` h --- Absorb a term from x-ingestX :: Bihom -> Maybe Integer -> Bihom-ingestX (a, b, c, d,-         e, f, g, h) (Just p)  = (b, a+b*p, d, c+d*p,-                                  f, e+f*p, h, g+h*p)-ingestX (_a, b, _c, d,-         _e, f, _g, h) Nothing = (b, b, d, d,-                                  f, f, h, h)--- Absorb a term from y-ingestY :: Bihom -> Maybe Integer -> Bihom-ingestY (a, b, c, d,-         e, f, g, h) (Just q)  = (c, d, a+c*q, b+d*q,-                                  g, h, e+g*q, f+h*q)-ingestY (_a, _b, c, d,-         _e, _f, g, h) Nothing = (c, d, c, d,-                                  g, h, g, h)---- Decide which of x and y to pull a term from-shouldIngestX :: Bihom -> Bool-shouldIngestX (_, _, _, _,-               0, 0, _, _) = False-shouldIngestX (_, _, _, _,-               0, _, 0, _) = True-shouldIngestX (a, b, c, _,-               e, f, g, _) = abs (g*e*b - g*a*f) > abs (f*e*c - g*a*f)+hom :: (Ord a, Num a, HasFractionField a, RealFrac (FractionField a)) => Hom a -> CF' a -> CF+hom (_n0, _n1,+     0,   0) _  = CF []+hom (_n0, _n1,+     0,   _d1) (CF []) = CF []+hom (n0, _n1,+     d0, _d1) (CF []) = valueToCF (insert n0 / insert d0)+hom h (CF (x:xs)) = case existsEmittable $ boundHom h (primitiveBound x) of+                     Just n ->  CF $ n : rest+                       where (CF rest) = hom (homEmit h (fromInteger n)) (CF (x:xs))+                     Nothing -> hom (homAbsorb h x) (CF xs) --- Apply a bihom to two continued fractions-bihom' :: Bihom -> [Integer] -> [Integer] -> [Integer]-bihom' (_, _, _, _,-        0, 0, 0, 0) _ _ = []-bihom' (0, 0, 0, 0,-        _, _, _, _) _ _ = [0]-bihom' bh x y = case biemit bh of-                Just (next, d) -> d : bihom' next x y-                Nothing -> if shouldIngestX bh then-                             bihom' (ingestX bh (safeHead x)) (safeRest x) y-                           else-                             bihom' (ingestY bh (safeHead y)) x (safeRest y)+bihomEmit :: Num a => Bihom a -> a -> Bihom a+bihomEmit (n0, n1, n2, n3,+           d0, d1, d2, d3) x = (d0,        d1,        d2,        d3,+                                n0 - d0*x, n1 - d1*x, n2 - d2*x, n3 - d3*x) -bihom :: Bihom -> CF -> CF -> CF-bihom bh (CF x) (CF y) = CF $ bihom' bh x y+bihomAbsorbX :: Num a => Bihom a -> a -> Bihom a+bihomAbsorbX (n0, n1, n2, n3,+              d0, d1, d2, d3) x = (n0*x + n1, n0, n2*x + n3, n2,+                                   d0*x + d1, d0, d2*x + d3, d2) -sqrt2 :: CF-sqrt2 = CF $ 1 : repeat 2+bihomAbsorbY :: Num a => Bihom a -> a -> Bihom a+bihomAbsorbY (n0, n1, n2, n3,+              d0, d1, d2, d3) y = (n0*y + n2, n1*y + n3, n0, n1,+                                   d0*y + d2, d1*y + d3, d0, d1) -exp1 :: CF-exp1 = CF (2 : concatMap triple [1..])-  where triple n = [1, 2 * n, 1]+bihomSubstituteX :: (Num a, HasFractionField a) => Bihom a -> Extended (FractionField a) -> Hom a+bihomSubstituteX (n0, n1, n2, n3,+                  d0, d1, d2, d3) (Finite x) = (n0*num + n1*den, n2*num + n3*den,+                                                d0*num + d1*den, d2*num + d3*den)+  where (num, den) = extract x+bihomSubstituteX (n0, _n1, n2, _n3,+                  d0, _d1, d2, _d3) Infinity = (n0, n2,+                                                d0, d2) -instance Eq CF where-  x == y = compare x y == EQ+bihomSubstituteY :: (Num a, HasFractionField a) => Bihom a -> Extended (FractionField a) -> Hom a+bihomSubstituteY (n0, n1, n2, n3,+                  d0, d1, d2, d3) (Finite y) = (n0*num + n2*den, n1*num + n3*den,+                                                d0*num + d2*den, d1*num + d3*den)+  where (num, den) = extract y+bihomSubstituteY (n0, n1, _n2, _n3,+                  d0, d1, _d2, _d3) Infinity = (n0, n1,+                                                d0, d1) -instance Ord CF where-  -- As [..., n, 1] represents the same number as [..., n+1]-  compare (CF [x]) (CF [y, 1]) = compare x (y+1)-  compare (CF [x, 1]) (CF [y]) = compare (x+1) y-  compare (CF [x]) (CF [y]) = compare x y+boundBihom :: (Ord a, Num a, HasFractionField a, Eq (FractionField a), Ord (FractionField a)) =>+              Bihom a -> Interval (FractionField a) -> Interval (FractionField a) -> Interval (FractionField a)+boundBihom bh x@(Interval ix sx) y@(Interval iy sy) = r1 `mergeInterval` r2 `mergeInterval` r3 `mergeInterval` r4+  where r1 = boundHom (bihomSubstituteX bh ix) y+        r2 = boundHom (bihomSubstituteY bh iy) x+        r3 = boundHom (bihomSubstituteX bh sx) y+        r4 = boundHom (bihomSubstituteY bh sy) x -  compare (CF (x:_)) (CF [y]) = if x < y then LT else GT-  compare (CF [x]) (CF (y:_)) = if x > y then GT else LT+select :: (Ord a, Num a, HasFractionField a, Eq (FractionField a), Ord (FractionField a)) =>+          Bihom a -> Interval (FractionField a) -> Interval (FractionField a) -> Bool+select bh x@(Interval ix sx) y@(Interval iy sy) = intX `smallerThan` intY+  where intX = if r1 `smallerThan` r2 then r2 else r1+        intY = if r3 `smallerThan` r4 then r3 else r4+        r1 = boundHom (bihomSubstituteX bh ix) y+        r2 = boundHom (bihomSubstituteX bh sx) y+        r3 = boundHom (bihomSubstituteY bh iy) x+        r4 = boundHom (bihomSubstituteY bh sy) x -  compare (CF (x:xs)) (CF (y:ys)) = case compare x y of-                                     EQ -> opposite $ compare (CF xs) (CF ys)-                                     o  -> o-    where opposite LT = GT-          opposite EQ = EQ-          opposite GT = LT+bihom :: (Ord a, Num a, HasFractionField a, RealFrac (FractionField a))+         => Bihom a -> CF' a -> CF' a -> CF+bihom bh (CF []) y = hom (bihomSubstituteX bh Infinity) y+bihom bh x (CF []) = hom (bihomSubstituteY bh Infinity) x+bihom bh (CF (x:xs)) (CF (y:ys)) = case existsEmittable $ boundBihom bh (primitiveBound x) (primitiveBound y) of+                   Just n -> CF $ n : rest+                     where (CF rest) = bihom (bihomEmit bh (fromInteger n)) (CF (x:xs)) (CF (y:ys))+                   Nothing -> if select bh (primitiveBound x) (primitiveBound y) then+                                let bh' = bihomAbsorbX bh x in bihom bh' (CF xs) (CF (y:ys))+                              else+                                let bh' = bihomAbsorbY bh y in bihom bh' (CF (x:xs)) (CF ys)  instance Num CF where   (+) = bihom (0, 1, 1, 0,-               1, 0, 0, 0)-  (*) = bihom (0, 0, 0, 1,-               1, 0, 0, 0)-  (-) = bihom (0, 1, -1, 0,-               1, 0,  0, 0)+               0, 0, 0, 1)+  (-) = bihom (0, -1, 1, 0,+               0,  0, 0, 1)+  (*) = bihom (1, 0, 0, 0,+               0, 0, 0, 1) -  fromInteger i = CF [i]-  abs x = if x > 0 then-             x-          else-            -x-  signum x | x < 0  = -1-           | x == 0 = 0-           | x > 0 = 1+  fromInteger n = CF [n] -instance Enum CF where-  toEnum = fromInteger . fromIntegral-  fromEnum = floor+  signum x = case 0 `compare` x of+              EQ -> 0+              LT -> 1+              GT -> -1 +  abs x | x < 0     = -x+        | otherwise = x+ instance Fractional CF where-  (/) = bihom (0, 1, 0, 0,-               0, 0, 1, 0)+  (/) = bihom (0, 0, 1, 0,+               0, 1, 0, 0) -  recip (CF [1]) = CF [1]-  recip (CF (0:xs)) = CF xs-  recip (CF xs) = CF (0:xs)+  fromRational = valueToCF -  fromRational r = fromInteger n / fromInteger d-    where n = numerator r-          d = denominator r+base :: Integer+base = 10 +rationalDigits :: Rational -> [Integer]+rationalDigits 0 = []+rationalDigits r = let d = num `quot` den in+                   d : rationalDigits (fromInteger base * (r - fromInteger d))+  where num = numerator r+        den = denominator r++digits :: CF -> [Integer]+digits = go (1, 0, 0, 1)+  where go (0, 0, _, _) _ = []+        go (p, _, q, _) (CF []) = rationalDigits (p % q)+        go h (CF (c:cs)) = case intervalDigit $ boundHom h (primitiveBound c) of+                            Nothing -> let h' = homAbsorb h c in go h' (CF cs)+                            Just d  -> d : go (homEmitDigit h d) (CF (c:cs))+        homEmitDigit (n0, n1,+                      d0, d1) d = (base * (n0 - d0*d), base * (n1 - d1*d),+                                   d0,                 d1)++cfString :: CF -> String+cfString (CF []) = "Infinity"+cfString cf | cf < 0 = '-' : cfString (-cf)+cfString cf = case digits cf of+               []     -> "0"+               [i]    -> show i+               (i:is) -> show i ++ "." ++ concatMap show is++instance Show CF where+  show = take 50 . cfString++instance Eq CF where+  a == b = a `compare` b == EQ++instance Ord CF where+  a `compare` b = head $ catMaybes $ zipWith comparePosition (nthPrimitiveBounds a) (nthPrimitiveBounds b)+ instance Real CF where-  -- Just take a pretty good rational approximation-  toRational cf = last $ take 20 (convergents cf)+  toRational = error "CF: toRational"  instance RealFrac CF where-  properFraction (CF [i]) = (fromIntegral i, 0)-  properFraction cf | cf < 0 = case properFraction (-cf) of-                                (b, a) -> (-b, -a)-  properFraction (CF (i:r)) = (fromIntegral i, CF r)+  properFraction cf = head $ mapMaybe checkValid $ nthPrimitiveBounds cf+    where checkValid (Interval (Finite a) (Finite b)) = if a <= b && truncate a == truncate b then+                                                          Just (truncate a, cf - fromInteger (truncate a))+                                                        else+                                                          Nothing+          checkValid _ = Nothing++cfcf :: CF' CF -> CF+cfcf = hom (1, 0, 0, 1)++instance Floating CF where+  exp r = cfcf (CF $ 1 : concatMap go [0..])+    where go n = [fromInteger (4*n+1) / r,+                  -2,+                  -fromInteger (4*n+3) / r,+                  2]++  -- TODO: restrict range+  log r = cfcf (CF $ 0 : concatMap go [0..])+    where go n = [fromInteger (2*n+1) / (r-1),+                  fromRational $ 2 % (n+1)]++  tan r = cfcf (CF $ 0 : concatMap go [0..])+    where go n = [fromInteger (4*n+1) / r,+                  -fromInteger (4*n+3) / r]++  sin r = bihom (0,2,0,0,+                 1,0,0,1) tanhalf tanhalf+    where tanhalf = tan (r / 2)++  cos r = bihom (-1,0,0,1,+                  1,0,0,1) tanhalf tanhalf+    where tanhalf = tan (r / 2)++  sinh r = bihom (1,0,0,-1,+                  0,1,1, 0) expr expr+    where expr = exp r++  cosh r = bihom (1,0,0,1,+                  0,1,1,0) expr expr+    where expr = exp r++  tanh r = bihom (1,0,0,-1,+                  1,0,0, 1) expr expr+    where expr = exp r

@@ -0,0 +1,163 @@+{-# LANGUAGE FlexibleInstances #-}+module Math.ContinuedFraction.Interval where++import Data.Ratio+import Numeric++data Extended a = Finite a | Infinity deriving (Eq)++data Interval a = Interval (Extended a) (Extended a) deriving (Eq)++instance Show (Interval Rational) where+  show (Interval a b) = "(" ++ showE a ++ ", " ++ showE b ++ ")"+    where showE Infinity = "Infinity"+          showE (Finite r) = show (fromRat r)++instance Num a => Num (Extended a) where+  Finite a + Finite b = Finite (a + b)+  Infinity + Finite _ = Infinity+  Finite _ + Infinity = Infinity+  Infinity + Infinity = error "Infinity + Infinity"++  Finite a * Finite b = Finite (a * b)+  Infinity * Finite a = Infinity+  -- Infinity * Finite a | a == 0 = error "Infinity * 0"+  --                     | otherwise = Infinity+  Finite a * i = i * Finite a+  Infinity * Infinity = undefined "Infinity * Infinity"++  negate (Finite r) = Finite (-r)+  negate Infinity = Infinity++  signum (Finite r) = Finite $ signum r+  signum Infinity = error "signum Infinity"++  abs (Finite r) = Finite $ abs r+  abs Infinity = Infinity++  fromInteger = Finite . fromInteger++instance (Show a) => Show (Extended a) where+  show (Finite r) = show r+  show Infinity = "Infinity"++smallerThan :: (Num a, Ord a) => Interval a -> Interval a -> Bool+Interval _ _ `smallerThan` Interval Infinity Infinity = False -- TODO CHECK+Interval Infinity Infinity `smallerThan` Interval _ _ = True+Interval (Finite a) Infinity `smallerThan` Interval (Finite b) Infinity = a >= b+Interval (Finite a) Infinity `smallerThan` Interval Infinity (Finite b) = a >= -b+Interval Infinity (Finite a) `smallerThan` Interval (Finite b) Infinity = a <= -b+Interval Infinity (Finite a) `smallerThan` Interval Infinity (Finite b) = a <= b+Interval (Finite i1) (Finite s1) `smallerThan` Interval Infinity (Finite _) = i1 <= s1+Interval (Finite i1) (Finite s1) `smallerThan` Interval (Finite _) Infinity = i1 <= s1+Interval Infinity (Finite _) `smallerThan` Interval (Finite i2) (Finite s2) = i2 > s2+Interval (Finite _) Infinity `smallerThan` Interval (Finite i2) (Finite s2) = i2 > s2+-- TODO: cache some of these comparisons+Interval (Finite i1) (Finite s1) `smallerThan` Interval (Finite i2) (Finite s2)+  =    (i1 <= s1 && i2 <= s2 && s1 - i1 <= s2 - i2)+    || (i1 >  s1 && i2 >  s2 && i1 - s1 >= i2 - s2)+    || (i1 <= s1 && i2 >  s2)++epsilon :: Rational+epsilon = 1 % 10^10++comparePosition :: Interval Rational -> Interval Rational -> Maybe Ordering+Interval (Finite i1) (Finite s1) `comparePosition` Interval (Finite i2) (Finite s2)+  | i1 > s1 = Nothing+  | i2 > s2 = Nothing+  | s1 < i2 = Just LT+  | s2 < i1 = Just GT+  | (s1 - i1) < epsilon && (s2 - i2) < epsilon = Just EQ+_ `comparePosition` _ = Nothing++intervalDigit :: (RealFrac a) => Interval a -> Maybe Integer+intervalDigit (Interval (Finite i) (Finite s)) = if i <= s && floor i == floor s && floor i >= 0 then+                                                   Just $ floor i+                                                 else+                                                   Nothing+intervalDigit _ = Nothing++subset :: Ord a => Interval a -> Interval a -> Bool+Interval _ _ `subset` Interval Infinity Infinity = True+Interval Infinity Infinity `subset` Interval _ _ = False+Interval Infinity (Finite s1) `subset` Interval Infinity (Finite s2) = s1 <= s2+Interval (Finite i1) Infinity `subset` Interval (Finite i2) Infinity = i1 >= i2+Interval Infinity (Finite _) `subset` Interval (Finite _) Infinity = False+Interval (Finite _) Infinity `subset` Interval Infinity (Finite _) = False+Interval (Finite i1) (Finite s1) `subset` Interval Infinity (Finite s2)+  | i1 <= s1 && s1 <= s2 = True+  | otherwise            = False+Interval (Finite i1) (Finite s1) `subset` Interval (Finite i2) Infinity+  | i1 <= s1 && i2 <= i1 = True+  | otherwise            = False+Interval Infinity (Finite s1) `subset` Interval (Finite i2) (Finite s2)+  | i2 > s2 && s1 <= s2 = True+  | otherwise            = False+Interval (Finite i1) Infinity `subset` Interval (Finite i2) (Finite s2)+  | i2 > s2 && i2 <= i1 = True+  | otherwise            = False+Interval (Finite i1) (Finite s1) `subset` Interval (Finite i2) (Finite s2)+  | i1 <= s1 && i2 <= s2 &&+    i2 <= i1 && s1 <= s2     = True+  | s1 <  i1 && s2 <  i2 &&+    i2 <= i1 && s1 <= s2     = True+  | i1 <= s1 && s2 <  i2 &&+    i2 <= i1 && i2 <= s1     = True+  | i1 <= s1 && s2 <  i2 &&+    i1 <= s2 && s1 <= s2     = True+  | otherwise                = False++elementOf :: (Ord a) => Extended a -> Interval a -> Bool+Infinity `elementOf` (Interval Infinity Infinity) = True+(Finite _) `elementOf` (Interval Infinity Infinity) = True+Infinity `elementOf` (Interval (Finite _) Infinity) = True+(Finite x) `elementOf` (Interval (Finite a) Infinity) = x >= a+Infinity `elementOf` (Interval Infinity (Finite _)) = True+(Finite x) `elementOf` (Interval Infinity (Finite b)) = x <= b+Infinity `elementOf` (Interval (Finite i) (Finite s)) = i > s+(Finite x) `elementOf` (Interval (Finite i) (Finite s))+  | i <= s = i <= x && x <= s+  | i >  s = i <= x || x <= s+  | otherwise = error "The impossible happened in elementOf"++-- Here we interpret Interval Infinity Infinity as the whole real line+mergeInterval :: (Ord a) => Interval a -> Interval a -> Interval a+mergeInterval (Interval Infinity Infinity) (Interval Infinity Infinity)+  = Interval Infinity Infinity+mergeInterval (Interval (Finite i) Infinity) (Interval Infinity Infinity)+  = Interval Infinity Infinity+mergeInterval (Interval Infinity (Finite s)) (Interval Infinity Infinity)+  = Interval Infinity Infinity+mergeInterval (Interval (Finite i) (Finite s)) (Interval Infinity Infinity)+  = Interval Infinity Infinity+mergeInterval (Interval Infinity (Finite s)) (Interval (Finite i) Infinity)+  | s >= i    = Interval Infinity Infinity+  | otherwise = Interval (Finite i) (Finite s)+mergeInterval (Interval Infinity (Finite s1)) (Interval Infinity (Finite s2))+  = Interval Infinity (Finite $ max s1 s2)+mergeInterval (Interval (Finite i1) Infinity) (Interval (Finite i2) Infinity)+  = Interval Infinity (Finite $ min i1 i2)+mergeInterval (Interval (Finite i1) (Finite s1)) (Interval (Finite i2) Infinity)+  | i1 <= s1 = Interval (Finite $ min i1 i2) Infinity+  | i1 >  s1 && i1 <= i2 = Interval (Finite i1) (Finite s1)+  | i1 >  s1 && i2 <= s1 = Interval Infinity Infinity+  | i1 >  s1 && i2 >  s1 = Interval (Finite i2) (Finite s1)+mergeInterval (Interval (Finite i1) (Finite s1)) (Interval Infinity (Finite s2))+  | i1 <= s1 = Interval Infinity (Finite $ max s1 s2)+  | i1 >  s1 && s2 <= s1 = Interval (Finite i1) (Finite s1)+  | i1 >  s1 && i1 <= s2 = Interval Infinity Infinity+  | i1 >  s1 && i1 >  s2 = Interval (Finite i1) (Finite s2)+mergeInterval int1@(Interval (Finite i1) (Finite s1)) int2@(Interval (Finite i2) (Finite s2))+  | i1 <= s1 && i2 <= s2 = Interval (Finite $ min i1 i2) (Finite $ max s1 s2)+  | i1 >  s1 && i2 >  s2 && (i1 <= s2 || i2 <= s1) = Interval Infinity Infinity+  | i1 >  s1 && i2 >  s2 = Interval (Finite $ min i1 i2) (Finite $ max s1 s2)+  | i1 >  s1 && i2 <= s2 = doTricky int2 int1+  | i1 <= s1 && i2 >  s2 = doTricky int1 int2+  | otherwise = error "The impossible happened in mergeInterval"+  where doTricky int1@(Interval (Finite i1) (Finite s1)) int2@(Interval (Finite i2) (Finite s2))+          | int1 `subset` int2         = int2+          | i2 <= s1 && i1 <= s2       = Interval Infinity Infinity+          | s1 < i2  = Interval (Finite i2) (Finite s1)+          | s2 < i1  = Interval (Finite i1) (Finite s2)+          | otherwise = error "The impossible happened in mergeInterval"+mergeInterval int1 int2 = mergeInterval int2 int1

@@ -0,0 +1,201 @@+module Math.ContinuedFraction.Simple+  (+    CF,+    digits,+    showCF,+    sqrt2,+    exp1+  ) where++import Data.Ratio++newtype CF = CF [Integer]++-- The coefficients of the homographic function (ax + b) / (cx + d)+type Hom = (Integer, Integer,+            Integer, Integer)++-- Possibly output a term+homEmittable :: Hom -> Maybe Integer+homEmittable (a, b,+              c, d) = if c /= 0 && d /= 0 && r == s then+                        Just r+                      else+                        Nothing+  where r = a `quot` c+        s = b `quot` d++homEmit :: Hom -> Integer -> Hom+homEmit (n0, n1,+         d0, d1) x = (d0,        d1,+                      n0 - d0*x, n1 - d1*x)++homAbsorb :: Hom -> Integer -> Hom+homAbsorb (n0, n1,+           d0, d1) x = (n0*x + n1, n0,+                        d0*x + d1, d0)++-- Apply a hom to a continued fraction+hom :: Hom -> CF -> CF+hom (0, 0,+     _, _) _ = CF [0]+hom (_, _,+     0, 0) _ = CF []+hom (n0, _,+     d0, _) (CF []) = fromRational (n0 % d0)+hom h (CF (x:xs)) = case homEmittable h of+                     Just d -> let (CF rest) = hom (homEmit h d) (CF (x:xs)) in CF (d : rest)+                     Nothing -> hom (homAbsorb h x) (CF xs)++-- The coefficients of the bihomographic function (axy + by + cx + d) / (exy + fy + gx + h)+type Bihom = (Integer, Integer, Integer, Integer,+              Integer, Integer, Integer, Integer)++bihomEmittable :: Bihom -> Maybe Integer+bihomEmittable (a, b, c, d,+                e, f, g, h) = if e /= 0 && f /= 0 && g /= 0 && h /= 0 && ratiosAgree then+                                Just r+                              else+                                Nothing+  where r = a `quot` e+        ratiosAgree = r == b `quot` f && r == c `quot` g && r == d `quot` h++bihomEmit :: Bihom -> Integer -> Bihom+bihomEmit (n0, n1, n2, n3,+           d0, d1, d2, d3) x = (d0,        d1,        d2,        d3,+                                n0 - d0*x, n1 - d1*x, n2 - d2*x, n3 - d3*x)++bihomAbsorbX :: Bihom -> Integer -> Bihom+bihomAbsorbX (n0, n1, n2, n3,+              d0, d1, d2, d3) x = (n0*x + n1, n0, n2*x + n3, n2,+                                   d0*x + d1, d0, d2*x + d3, d2)++bihomAbsorbY :: Bihom -> Integer -> Bihom+bihomAbsorbY (n0, n1, n2, n3,+              d0, d1, d2, d3) y = (n0*y + n2, n1*y + n3, n0, n1,+                                   d0*y + d2, d1*y + d3, d0, d1)++-- Decide which of x and y to pull a term from+shouldIngestX :: Bihom -> Bool+shouldIngestX (_, _, _, _,+               _, 0, _, 0) = True+shouldIngestX (_, _, _, _,+               _, _, 0, 0) = False+shouldIngestX (_a, b, c, d,+               _e, f, g, h) = abs (g*h*b - g*d*f) < abs (f*h*c - g*d*f)++-- Apply a bihom to two continued fractions+bihom :: Bihom -> CF -> CF -> CF+bihom (_, _, _, _,+       0, 0, 0, 0) _ _ = CF []+bihom (0, 0, 0, 0,+       _, _, _, _) _ _ = CF [0]+bihom (n0, _n1, n2, _n3,+       d0, _d1, d2, _d3) (CF []) y = hom (n0, n2,+                                          d0, d2) y+bihom (n0, n1, _n2, _n3,+       d0, d1, _d2, _d3) x (CF []) = hom (n0, n1,+                                          d0, d1) x+bihom bh (CF (x:xs)) (CF (y:ys)) = case bihomEmittable bh of+                                    Just d -> CF $ d : rest+                                      where (CF rest) = bihom (bihomEmit bh d) (CF (x:xs)) (CF (y:ys))+                                    Nothing -> if shouldIngestX bh then+                                                 bihom (bihomAbsorbX bh x) (CF xs) (CF (y:ys))+                                               else+                                                 bihom (bihomAbsorbY bh y) (CF (x:xs)) (CF ys)++sqrt2 :: CF+sqrt2 = CF $ 1 : repeat 2++exp1 :: CF+exp1 = CF (2 : concatMap triple [1..])+  where triple n = [1, 2 * n, 1]++instance Eq CF where+  x == y = compare x y == EQ++instance Ord CF where+  -- As [..., n, 1] represents the same number as [..., n+1]+  compare (CF [x]) (CF [y, 1]) = compare x (y+1)+  compare (CF [x, 1]) (CF [y]) = compare (x+1) y+  compare (CF [x]) (CF [y]) = compare x y++  compare (CF (x:_)) (CF [y]) = if x < y then LT else GT+  compare (CF [x]) (CF (y:_)) = if x > y then GT else LT++  compare (CF (x:xs)) (CF (y:ys)) = case compare x y of+                                     EQ -> opposite $ compare (CF xs) (CF ys)+                                     o  -> o+    where opposite LT = GT+          opposite EQ = EQ+          opposite GT = LT++instance Num CF where+  (+) = bihom (0, 1, 1, 0,+               0, 0, 0, 1)+  (*) = bihom (1, 0, 0, 0,+               0, 0, 0, 1)+  (-) = bihom (0, -1, 1, 0,+               0,  0, 0, 1)++  fromInteger i = CF [i]+  abs x = if x > 0 then+             x+          else+            -x+  signum x | x < 0  = -1+           | x == 0 = 0+           | x > 0 = 1+++instance Fractional CF where+  (/) = bihom (0, 0, 1, 0,+               0, 1, 0, 0)++  recip (CF [1]) = CF [1]+  recip (CF (0:xs)) = CF xs+  recip (CF xs) = CF (0:xs)++  fromRational r = if rest == 0 then+                CF [d]+              else+                let (CF ds)  = fromRational (recip rest) in CF (d:ds)+    where (d, rest) = properFraction r++instance Real CF where+  toRational _ = undefined++instance RealFrac CF where+  properFraction (CF [i]) = (fromIntegral i, 0)+  properFraction cf | cf < 0 = case properFraction (-cf) of+                                (b, a) -> (-b, -a)+  properFraction (CF (i:r)) = (fromIntegral i, CF r)++rationalDigits :: Rational -> [Integer]+rationalDigits 0 = []+rationalDigits r = let d = num `quot` den in+                   d : rationalDigits (fromInteger 10 * (r - fromInteger d))+  where num = numerator r+        den = denominator r++digits :: CF -> [Integer]+digits = go (1, 0, 0, 1)+  where go (0, 0, _, _) _ = []+        go (p, _, q, _) (CF []) = rationalDigits (p % q)+        go h (CF (c:cs)) = case homEmittable h of+                            Nothing -> let h' = homAbsorb h c in go h' (CF cs)+                            Just d  -> d : go (homEmitDigit h d) (CF (c:cs))+        homEmitDigit (n0, n1,+                      d0, d1) d = (10 * (n0 - d0*d), 10 * (n1 - d1*d),+                                   d0,               d1)++-- | Produce a decimal representation of a number+showCF :: CF -> String+showCF cf | cf < 0 = "-" ++ show (-cf)+showCF (CF [i])   = show i+showCF (CF (i:r)) = show i ++ "." ++ decimalDigits+  where decimalDigits = concatMap show $ tail $ digits (CF (0:r))++-- Should make this cleverer+instance Show CF where+  show = take 15 . showCF

@@ -0,0 +1,48 @@+{-# LANGUAGE TemplateHaskell, FlexibleInstances #-}+module Main where++import Data.Maybe+import Data.Ratio++import Math.ContinuedFraction+import Math.ContinuedFraction.Interval++import Test.QuickCheck+import Test.QuickCheck.Function+import Test.Framework.TH+import Test.Framework.Providers.QuickCheck2++instance Arbitrary (Extended Rational) where+  arbitrary = do+    b <- arbitrary :: Gen Bool+    if b then+      return Infinity+    else do+      n <- choose (-10, 10)+      return $ Finite (n % 1)++instance Arbitrary (Interval Rational) where+  arbitrary = do+    (i, s) <- suchThat arbitrary (\(i,s) -> i /= s) :: Gen (Extended Rational, Extended Rational)+    return $ Interval i s++prop_sensiblePrimitiveBound x = fromInteger x `elementOf` primitiveBound x+  where types = x :: Integer++prop_sensibleMergeInterval a b = a `subset` ab && b `subset` ab+  where types = (a :: Interval Rational, b :: Interval Rational)+        ab = a `mergeInterval` b++finitePrimitiveBounds (CF cf) = zipWith boundHom homs (map primitiveBound cf)+  where homs = scanl homAbsorb (1,0,0,1) cf++prop_primitiveBoundsContain a b = all ((Finite $ a + b) `elementOf`) $ finitePrimitiveBounds (valueToCF a + valueToCF b)+  where types = (a :: Rational, b :: Rational)++prop_sensibleEuclidean i = case existsEmittable i of+                            Just n -> i `subset` primitiveBound n+                            Nothing -> True+  where types = i :: Interval Rational++main :: IO ()+main = $defaultMainGenerator