{-# LANGUAGE MagicHash #-}
{-# LANGUAGE MultiWayIf #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE DataKinds #-}
module Data.CReal.Internal
( CReal(..)
, atPrecision
, crealPrecision
, shiftL
, shiftR
, (/.)
, log2
, log10
, decimalDigitsAtPrecision
, rationalToDecimal
) where
import Data.Ratio (numerator,denominator,(%))
import GHC.Base (Int(..))
import GHC.Integer.Logarithms (integerLog2#, integerLogBase#)
import GHC.TypeLits
infixl 7 /.
default ()
newtype CReal (n :: Nat) = CR (Int -> Integer)
crealPrecision :: KnownNat n => CReal n -> Int
crealPrecision = fromInteger . natVal
atPrecision :: CReal n -> Int -> Integer
(CR x) `atPrecision` p = x p
instance KnownNat n => Show (CReal n) where
show x = showAtPrecision (crealPrecision x) x
instance Num (CReal n) where
fromInteger i = CR (\p -> i * 2 ^ p)
negate (CR x) = CR (negate . x)
abs (CR x) = CR (abs . x)
{-# INLINE (+) #-}
CR x1 + CR x2 = CR (\p -> let n1 = x1 (p + 2)
n2 = x2 (p + 2)
in (n1 + n2) /. 4)
{-# INLINE (*) #-}
CR x1 * CR x2 = CR (\p -> let s1 = log2 (abs (x1 0) + 2) + 3
s2 = log2 (abs (x2 0) + 2) + 3
n1 = x1 (p + s2)
n2 = x2 (p + s1)
in (n1 * n2) /. 2^(p + s1 + s2) )
signum x = CR (\p -> signum (x `atPrecision` p) * 2^p)
instance Fractional (CReal n) where
fromRational n = fromInteger (numerator n) / fromInteger (denominator n)
{-# INLINE recip #-}
recip (CR x) = CR (\p -> let s = findFirstMonotonic ((3 <=) . abs . x)
n = x (p + 2 * s + 2)
in 2^(2 * p + 2 * s + 2) /. n)
instance KnownNat n => Eq (CReal n) where
x == y = let p = crealPrecision x
in (x - y) `atPrecision` p == 0
instance KnownNat n => Ord (CReal n) where
compare x y = let p = crealPrecision x
in compare ((x - y) `atPrecision` p) 0
max (CR x) (CR y) = CR (\p -> max (x p) (y p))
min (CR x) (CR y) = CR (\p -> min (x p) (y p))
shiftR :: CReal n -> Int -> CReal n
shiftR (CR x) n = CR (\p -> let p' = p - n
in if p' >= 0
then x p'
else x 0 /. 2^(-p'))
shiftL :: CReal n -> Int -> CReal n
shiftL x = shiftR x . negate
showAtPrecision :: Int -> CReal n -> String
showAtPrecision p (CR x) = let places = decimalDigitsAtPrecision p
r = x p % 2^p
in rationalToDecimal places r
decimalDigitsAtPrecision :: Int -> Int
decimalDigitsAtPrecision 0 = 0
decimalDigitsAtPrecision p = log10 (2^p) + 1
rationalToDecimal :: Int -> Rational -> String
rationalToDecimal places r = s
where ds = show ((numerator r * 10^places) /. denominator r)
(is, fs) = splitAt (length ds - places) ds
is' = case is of
"" -> "0"
"-" -> "-0"
_ -> is
suff = case places of
0 -> ""
_ -> '.' : take places (fs ++ repeat '0')
s = is' ++ suff
(/.) :: Integer -> Integer -> Integer
n /. d = round (n % d)
log2 :: Integer -> Int
log2 x = I# (integerLog2# x)
log10 :: Integer -> Int
log10 x = I# (integerLogBase# 10 x)
findFirstMonotonic :: (Int -> Bool) -> Int
findFirstMonotonic p = binarySearch l' u'
where (l', u') = findBounds 0 1
findBounds l u = if p u then (l, u)
else findBounds u (u*2)
binarySearch l u = let m = l + ((u - l) `div` 2)
in if | l+1 == u -> l
| p m -> binarySearch l m
| otherwise -> binarySearch m u