From d57ae9b4a00eaf0aa02ce463f73d2645f3ba0d3d Mon Sep 17 00:00:00 2001 From: rnhmjoj Date: Thu, 20 Aug 2015 20:40:13 +0200 Subject: [PATCH] Rewrite arithmetics --- src/Data/Number/Instances.hs | 19 +---- src/Data/Number/Internal.hs | 145 ++++++++++++++++++++++++++--------- 2 files changed, 113 insertions(+), 51 deletions(-) diff --git a/src/Data/Number/Instances.hs b/src/Data/Number/Instances.hs index 1ff0f53..39e036f 100644 --- a/src/Data/Number/Instances.hs +++ b/src/Data/Number/Instances.hs @@ -28,9 +28,9 @@ instance Traversable Continued where -- | The sign is given by the first number of the fraction. -- Other number are always considered positive. instance Num Number where - (+) = operator (0, 1, 1, 0, 1, 0, 0, 0) - (-) = operator (0, 1, -1, 0, 1, 0, 0, 0) - (*) = operator (0, 0, 0, 1, 1, 0, 0, 0) + (+) = biHom (0, 1, 1, 0, 1, 0, 0, 0) + (-) = biHom (0, 1, -1, 0, 1, 0, 0, 0) + (*) = biHom (0, 0, 0, 1, 1, 0, 0, 0) abs (M x) = x abs x = x @@ -50,7 +50,7 @@ instance Real Number where -- | Allows division between 'Number's and conversion from a rational instance Fractional Number where - (/) = operator (0, 1, 0, 0, 0, 0, 1, 0) + (/) = biHom (0, 1, 0, 0, 0, 0, 1, 0) fromRational = toNumber @@ -62,14 +62,3 @@ fromNumber E = 0 fromNumber (M x) = negate (fromNumber x) fromNumber (x :| E) = fromIntegral x fromNumber (x :| xs) = fromIntegral x + 1 / (fromNumber xs) - --- | Convert a 'RealFrac' number into a 'Number' -toNumber :: (Show a, RealFrac a) => a -> Number -toNumber 0 = E -toNumber x0 - | x0 < 0 = M (n :| toNumber x1) - | otherwise = n :| toNumber x1 - where - (n,f) = properFraction (abs x0) - x1 | f < 1e-6 = 0 - | otherwise = 1/f diff --git a/src/Data/Number/Internal.hs b/src/Data/Number/Internal.hs index 87abc4a..f36534a 100644 --- a/src/Data/Number/Internal.hs +++ b/src/Data/Number/Internal.hs @@ -1,10 +1,12 @@ -- | Data.Number internals module Data.Number.Internal -( Matrix -, operator +( Hom, BiHom +, hom, biHom +, toNumber , cut , first , rest +, join , split ) where @@ -12,50 +14,114 @@ import Data.Number.Types import Data.Number.Peano import Data.Ratio +type Hom = (Whole, Whole, Whole, Whole) +type BiHom = (Whole, Whole, Whole, Whole, + Whole, Whole, Whole, Whole) --- | Operator Matrix -type Matrix = (Whole, Whole, Whole, Whole, Whole, Whole, Whole, Whole) - --- | Continued fraction operator (implements Gosper's arithmetics) +-- | Homographic function -- --- Given two 'Number' @x@, @y@ and the operator matrix +-- Given the 'Hom' matrix -- --- <> +-- <> -- --- calculates +-- and a 'Number' @x@ calculates -- --- <> +-- <> -- -- See for a complete -- explanation. -operator :: Matrix -> Number -> Number -> Number -operator c x y = - case operator' c x y False of - [] -> E - m -> if head m < 0 - then M $ fromList (map toNat m) - else fromList (map toNat m) - where - fromList [] = E - fromList (x:xs) = x :| fromList xs +hom :: Hom -> Number -> Number +hom (0, 0, _, _) _ = E +hom (a, _, c, _) E = toNumber (fromPeano a % fromPeano c) +hom h x = case maybeEmit h of + Just d -> join d (hom (emit h d) x) + Nothing -> hom (absorb h x0) x' + where (x0, x') = split x -operator' :: Matrix -> Number -> Number -> Bool -> [Whole] -operator' (_,_,_,_,0,0,0,0) _ _ _ = [] -operator' (a,b,c,d,e,f,g,h) x y end - | t = r : operator' (e, f, g, h, a-e*r, b-f*r, c-g*r, d-h*r) x y end - | x/=E && s = operator' (b, a+b*p, d, c+d*p, f, e+f*p, h, g+h*p) x' y end - | x==E && s = operator' (b, b, d, d, f, f, h, h) E y end - | y/=E = operator' (c, d, a+c*q, b+d*q, g, h, e+g*q, f+h*q) x y' end - | otherwise = operator' (c, d, c, d, g, h, g, h) x E True - where - r = a // e - (p, x') = split x - (q, y') = split y - t = not (any (==0) [e,f,g,h]) && all (==r) [b//f, c//g, d//h] - s | end = True - | any (==0) [f,g,e,h] = False - | otherwise = abs (b%f - a%e) > abs (c%g - a%e) +-- Homographic helpers -- + +maybeEmit :: Hom -> Maybe Whole +maybeEmit (a, b, c, d) = + if c /= 0 && d /= 0 && r == s + then Just r + else Nothing + where r = a // c + s = b // d + + +emit :: Hom -> Whole -> Hom +emit (a, b, c, d) x = (c, d, a - c*x, b - d*x) + + +absorb :: Hom -> Whole -> Hom +absorb (a, b, c, d) x = (a*x + b, a, c*x + d, c) + + +-- | Bihomographic function +-- +-- Given a 'Hom' matrix +-- +-- <> +-- +-- and two 'Number' @x@ and @y@ calculates +-- +-- <> +biHom :: BiHom -> Number -> Number -> Number +biHom (0, 0, 0, 0, _, _, _, _) _ _ = E +biHom (a, _, c, _, e, _, g, _) E y = hom (a, c, e, g) y +biHom (a, b, _, _, e, f, _, _) x E = hom (a, b, e, f) x +biHom h x y = case maybeBiEmit h of + Just d -> join d (biHom (biEmit h d) x y) + Nothing -> if fromX h + then biHom (biAbsorbX h x0) x' y + else biHom (biAbsorbY h y0) x y' + where + (x0, x') = split x + (y0, y') = split y + + +-- Bihomographic helpers + +maybeBiEmit :: BiHom -> Maybe Whole +maybeBiEmit (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 = quot a e + ratiosAgree = r == b // f && r == c // g && r == d // h + +biEmit :: BiHom -> Whole -> BiHom +biEmit (a, b, c, d, + e, f, g, h) x = (e, f, g, h, + a - e*x, b - f*x, c - g*x, d - h*x) + +biAbsorbX :: BiHom -> Whole -> BiHom +biAbsorbX (a, b, c, d, + e, f, g, h) x = (a*x + b, a, c*x + d, c, + e*x + f, e, g*x + h, g) + +biAbsorbY :: BiHom -> Whole -> BiHom +biAbsorbY (a, b, c, d, + e, f, g, h) y = (a*y + c, b*y + d, a, b, + e*y + g, f*y + h, e, f) + + +fromX :: BiHom -> Bool +fromX (_, _, _, _, _, 0, _, 0) = True +fromX (_, _, _, _, _, _, 0, 0) = False +fromX (_, b, c, d, _, f, g, h) = abs (g*h*b - g*d*f) < abs (f*h*c - g*d*f) + + +-- | Convert a 'RealFrac' number into a 'Number' +toNumber :: RealFrac a => a -> Number +toNumber 0 = E +toNumber x + | x < 0 = M (toNumber (-x)) + | x' == 0 = x0 :| E + | otherwise = x0 :| toNumber (recip x') + where (x0, x') = properFraction x -- | Truncate a 'Number' to a given length @n@ cut :: Nat -> Number -> Number @@ -72,6 +138,12 @@ split :: Number -> (Whole, Number) split x = (first x, rest x) +-- | Essentially the inverse of split +join :: Whole -> Number -> Number +join (Whole x0 Neg) = M . (x0 :|) +join (Whole x0 Pos) = (x0 :|) + + -- | Extract the first natural of the fraction as a 'Whole' number first :: Number -> Whole first E = 0 @@ -81,6 +153,7 @@ first (x:|_) = Whole x Pos -- | Extract the "tail" of a 'Number' as a new 'Number' +-- -- Equivalent to @(x - floor x)@ for a floating point. rest :: Number -> Number rest E = E