diff --git a/cf.cabal b/cf.cabal --- a/cf.cabal +++ b/cf.cabal @@ -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 diff --git a/src/Math/ContinuedFraction.hs b/src/Math/ContinuedFraction.hs --- a/src/Math/ContinuedFraction.hs +++ b/src/Math/ContinuedFraction.hs @@ -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 diff --git a/src/Math/ContinuedFraction/Interval.hs b/src/Math/ContinuedFraction/Interval.hs new file mode 100644 --- /dev/null +++ b/src/Math/ContinuedFraction/Interval.hs @@ -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 diff --git a/src/Math/ContinuedFraction/Simple.hs b/src/Math/ContinuedFraction/Simple.hs new file mode 100644 --- /dev/null +++ b/src/Math/ContinuedFraction/Simple.hs @@ -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 diff --git a/tests/Tests.hs b/tests/Tests.hs new file mode 100644 --- /dev/null +++ b/tests/Tests.hs @@ -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