Rewrite arithmetics

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)
+  (+) = biHom (0, 1,  1, 0, 1, 0, 0, 0)
-  (-) = operator (0, 1, -1, 0, 1, 0, 0, 0)
+  (-) = biHom (0, 1, -1, 0, 1, 0, 0, 0)
-  (*) = operator (0, 0,  0, 1, 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

src/Data/Number/Internal.hs

@@ -1,10 +1,12 @@
 -- | Data.Number internals
 module Data.Number.Internal
-( Matrix
+( Hom, BiHom
-, operator
+, 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
+-- | Homographic function
-type Matrix = (Whole, Whole, Whole, Whole, Whole, Whole, Whole, Whole)
-
--- | Continued fraction operator (implements Gosper's arithmetics)
 --
--- Given two 'Number' @x@, @y@ and the operator matrix
+-- Given the 'Hom' matrix
 --
--- <<https://i.imgur.com/Hm7TiIH.png>>
+-- <<https://i.imgur.com/iGobkbj.png>>
 --
--- calculates 
+-- and a 'Number' @x@ calculates
 --
--- <<https://i.imgur.com/IZvQmy9.png>>
+-- <<https://i.imgur.com/pCq29U3.png>>
 --
 -- See <http://perl.plover.com/yak/cftalk/INFO/gosper.txt> for a complete
 -- explanation.
-operator :: Matrix -> Number -> Number -> Number
+hom :: Hom -> Number -> Number
-operator c x y =
+hom (0, 0, _, _) _ = E
-  case operator' c x y False of
+hom (a, _, c, _) E = toNumber (fromPeano a % fromPeano c)
-    [] -> E
+hom h x = case maybeEmit h of
-    m  -> if head m < 0
+  Just d  -> join d (hom (emit h d) x)
-      then M $ fromList (map toNat m)
+  Nothing -> hom (absorb h x0) x'
-      else     fromList (map toNat m)
+  where (x0, x') = split x
-  where
-    fromList []     = E
-    fromList (x:xs) = x :| fromList xs
 
-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
+--
+-- <<https://i.imgur.com/Hm7TiIH.png>>
+--
+-- and two 'Number' @x@ and @y@ calculates
+--
+-- <<https://i.imgur.com/IZvQmy9.png>>
+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