number/src/Data/Number/Internal.hs

-- | Data.Number internals
module Data.Number.Internal
( operator
, cut
, first
, rest
, split
) where

import Data.Number.Types
import Data.Number.Peano
import Data.Ratio


-- | Operator Matrix
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
-- @\<a, b, c, d, e, f, g, h\>@
-- calculates @z = (a + bx + cy + dxy) / (e + fx + gy + hxy)@
--
-- See <http://perl.plover.com/yak/cftalk/INFO/gosper.txt> 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

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)


-- | 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)


-- | 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