166 lines
4.0 KiB
Haskell
166 lines
4.0 KiB
Haskell
-- | Data.Number internals
|
|
module Data.Number.Internal
|
|
( Hom, BiHom
|
|
, hom, biHom
|
|
, toNumber
|
|
, cut
|
|
, first
|
|
, rest
|
|
, join
|
|
, split
|
|
) where
|
|
|
|
import Data.Number.Types
|
|
import Data.Number.Peano
|
|
import Data.Ratio
|
|
|
|
-- | Homographic function coefficients matrix
|
|
type Hom = (Whole, Whole, Whole, Whole)
|
|
|
|
-- | Bihomographic function coefficients matrix
|
|
type BiHom = (Whole, Whole, Whole, Whole,
|
|
Whole, Whole, Whole, Whole)
|
|
|
|
-- | Homographic function
|
|
--
|
|
-- Given the 'Hom' matrix
|
|
--
|
|
-- <<https://i.imgur.com/iGobkbj.png>>
|
|
--
|
|
-- and a 'Number' @x@ calculates
|
|
--
|
|
-- <<https://i.imgur.com/pCq29U3.png>>
|
|
--
|
|
-- See <http://perl.plover.com/yak/cftalk/INFO/gosper.txt> for a complete
|
|
-- explanation.
|
|
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
|
|
|
|
|
|
-- 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 'BiHom' 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
|
|
cut _ E = E
|
|
cut n (M x) = M (cut n x)
|
|
cut n _ | n <= 0 = E
|
|
cut n (x :| xs) = x :| cut (n-1) xs
|
|
|
|
|
|
-- | Split a Number into a 'Whole' (the most significant of the fraction)
|
|
-- and the rest of the Number. Equivalent to @(floor x, x - floor x)@
|
|
-- for a floating point.
|
|
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
|
|
first (M E) = 0
|
|
first (M (x:|_)) = Whole x Neg
|
|
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
|
|
rest (M E) = E
|
|
rest (M x) = M (rest x)
|
|
rest (_:|xs) = xs
|