module Numeric.Decimal.Number
( Sign(..)
, negateSign
, xorSigns
, Coefficient
, numDigits
, Exponent
, Payload
, Decimal(..)
, BasicDecimal
, ExtendedDecimal
, GeneralDecimal
, zero
, oneHalf
, one
, negativeOne
, two
, ten
, infinity
, qNaN
, sNaN
, flipSign
, cast
, fromBool
, fromOrdering
, isPositive
, isNegative
, isFinite
, isZero
, isNormal
, isSubnormal
, adjustedExponent
, integralValue
) where
import Prelude hiding (exponent)
import Control.DeepSeq (NFData(..))
import Control.Monad (join)
import Data.Bits (Bits(..), FiniteBits(..))
import Data.Char (isSpace)
import Data.Coerce (coerce)
import Data.Ratio (numerator, denominator, (%))
import Numeric.Natural (Natural)
import Text.ParserCombinators.ReadP (readP_to_S)
import {-# SOURCE #-} Numeric.Decimal.Arithmetic
import {-# SOURCE #-} Numeric.Decimal.Conversion
import Numeric.Decimal.Precision
import Numeric.Decimal.Rounding
import {-# SOURCE #-} qualified Numeric.Decimal.Operation as Op
import qualified GHC.Real
data Sign = Pos
| Neg
deriving (Eq, Enum)
instance NFData Sign where
rnf s = s `seq` ()
negateSign :: Sign -> Sign
negateSign Pos = Neg
negateSign Neg = Pos
xorSigns :: Sign -> Sign -> Sign
xorSigns Pos Pos = Pos
xorSigns Pos Neg = Neg
xorSigns Neg Pos = Neg
xorSigns Neg Neg = Pos
signFunc :: Num a => Sign -> a -> a
signFunc Pos = id
signFunc Neg = negate
signMatch :: (Num a, Eq a) => a -> Sign
signMatch x = case signum x of
-1 -> Neg
_ -> Pos
type Coefficient = Natural
type Exponent = Integer
type Payload = Coefficient
data Decimal p r
= Num { sign :: Sign
, coefficient :: Coefficient
, exponent :: Exponent
}
| Inf { sign :: Sign
}
| NaN { sign :: Sign
, signaling :: Bool
, payload :: Payload
}
type BasicDecimal = Decimal P9 RoundHalfUp
type ExtendedDecimal p = Decimal p RoundHalfEven
type GeneralDecimal = ExtendedDecimal PInfinite
instance Show (Decimal p r) where
showsPrec d n = showParen (d > 0 && isNegative n) $ toScientificString n
instance (Precision p, Rounding r) => Read (Decimal p r) where
readsPrec _ str = [ (cast n, s)
| (n, s) <- readParen False
(readP_to_S toNumber . dropWhile isSpace) str ]
decimalPrecision :: Decimal p r -> p
decimalPrecision = undefined
instance Precision p => Precision (Decimal p r) where
precision = precision . decimalPrecision
eMax = eMax . decimalPrecision
eMin = eMin . decimalPrecision
evalOp :: Arith p r a -> a
evalOp op = let Right r = evalArith op newContext in r
evalOp' :: Arith p RoundHalfEven a -> a
evalOp' = evalOp
compareDecimal :: Decimal a b -> Decimal c d -> Either GeneralDecimal Ordering
compareDecimal x y = evalOp (x `Op.compare` y)
instance Eq (Decimal p r) where
x == y = case compareDecimal x y of
Right EQ -> True
_ -> False
instance Ord (Decimal p r) where
compare x y = case compareDecimal x y of
Right o -> o
Left _ -> evalOp (x `Op.compareTotal` y)
x < y = case compareDecimal x y of
Right LT -> True
_ -> False
x <= y = case compareDecimal x y of
Right LT -> True
Right EQ -> True
_ -> False
x > y = case compareDecimal x y of
Right GT -> True
_ -> False
x >= y = case compareDecimal x y of
Right GT -> True
Right EQ -> True
_ -> False
max x y = evalOp (Op.max x y)
min x y = evalOp (Op.min x y)
instance (Precision p, Rounding r) => Enum (Decimal p r) where
succ x = evalOp (x `Op.add` one)
pred x = evalOp (x `Op.subtract` one)
toEnum = fromIntegral
fromEnum Num { sign = s, coefficient = c, exponent = e }
| e >= 0 = signFunc s (fromIntegral c * 10^ e )
| otherwise = signFunc s (fromIntegral $ c `quot` (10^(-e)))
fromEnum _ = 0
enumFrom x = enumFromWith x one
enumFromThen x y = let i = y - x
in x : enumFromWith y i
enumFromTo x z = takeWhile (<= z) $ enumFromWith x one
enumFromThenTo x y z = let i = y - x
cmp | i < 0 = (>=)
| otherwise = (<=)
in takeWhile (`cmp` z) $ x : enumFromWith y i
enumFromWith :: (Precision p, Rounding r)
=> Decimal p r -> Decimal p r -> [Decimal p r]
enumFromWith x i = x : enumFromWith (x + i) i
instance (Precision p, Rounding r) => Num (Decimal p r) where
x + y = evalOp (x `Op.add` y)
x - y = evalOp (x `Op.subtract` y)
x * y = evalOp (x `Op.multiply` y)
negate = evalOp . Op.minus
abs = evalOp . Op.abs
signum n = case n of
Num { sign = s, coefficient = c }
| c == 0 -> zero { sign = s }
| otherwise -> one { sign = s }
Inf { sign = s } -> one { sign = s }
_ -> qNaN
fromInteger x = cast
Num { sign = signMatch x
, coefficient = fromInteger (abs x)
, exponent = 0
}
instance (Precision p, Rounding r) => Real (Decimal p r) where
toRational Num { sign = s, coefficient = c, exponent = e }
| e >= 0 = fromInteger $ signFunc s (fromIntegral c * 10^e)
| otherwise = signFunc s (fromIntegral c) % 10^(-e)
toRational n = signFunc (sign n) $ case n of
Inf{} -> GHC.Real.infinity
_ -> GHC.Real.notANumber
instance (FinitePrecision p, Rounding r) => Fractional (Decimal p r) where
x / y = evalOp (x `Op.divide` y)
fromRational r = let n = fromInteger (numerator r) :: GeneralDecimal
d = fromInteger (denominator r) :: GeneralDecimal
in evalOp (n `Op.divide` d)
instance (FinitePrecision p, Rounding r) => RealFrac (Decimal p r) where
properFraction x@Num { sign = s, coefficient = c, exponent = e }
| e < 0 = (n, f)
| otherwise = (signFunc s (fromIntegral c * 10^e), zero)
where n = signFunc s (fromIntegral q)
f = x { coefficient = r }
(q, r) = c `quotRem` (10^(-e))
properFraction nan = (0, nan)
continuedFraction :: FinitePrecision p
=> Int -> Integer -> [(Integer, Integer)] -> ExtendedDecimal p
continuedFraction m b0 ((a1, b1) : ps) = convergent (max 0 $ m - 2) x0 x1 x1' ps
where x0 = (aa0, bb0)
x1 = (aa1, bb1)
x1' = fromRational (aa1 % bb1)
aa0 = b0
bb0 = 1
aa1 = b1 * b0 + a1
bb1 = b1
convergent m (aa0, bb0) x1@(aa1, bb1) x1' ((a2, b2) : ps)
| m == 0 && x2' == x1' = x2'
| otherwise = convergent (max 0 $ m - 1) x1 x2 x2' ps
where x2 = (aa2, bb2)
x2' = fromRational (aa2 % bb2)
aa2 = b2 * aa1 + a2 * aa0
bb2 = b2 * bb1 + a2 * bb0
convergent _ _ _ x [] = x
continuedFraction _ b0 [] = fromInteger b0
infiniteSeries :: FinitePrecision p
=> (ExtendedDecimal p -> ExtendedDecimal p -> ExtendedDecimal p)
-> [ExtendedDecimal p] -> ExtendedDecimal p
infiniteSeries op ~(x:xs) = series x xs
where series n (x:xs)
| n' == n = n'
| otherwise = series n' xs
where n' = n `op` x
series n [] = n
seriesArctan :: FinitePrecision p => Decimal p r -> ExtendedDecimal p
seriesArctan z = infiniteSeries (+) (z' : series True three z')
where series neg d z =
let z' = z * z2
n | neg = flipSign z'
| otherwise = z'
in (n / d) : series (not neg) (d + two) z'
z' = castRounding z
z2 = z' * z'
seriesArcsin :: FinitePrecision p => Decimal p r -> ExtendedDecimal p
seriesArcsin z = infiniteSeries (+) (z' : series one two z' three)
where series n d z i =
let z' = z * z2
in (n * z') / (d * i) : series (n * i) (d * (i + one)) z' (i + two)
z' = castRounding z
z2 = z' * z'
arctan :: (FinitePrecision p, Rounding r) => Decimal p r -> ExtendedDecimal p
arctan z@Num { } = arctan' (toRational z)
arctan Inf { sign = s } = signFunc s halfPi
arctan _ = qNaN
arctan' :: FinitePrecision p => Rational -> ExtendedDecimal p
arctan' z = continuedFraction m 0 $ (x, y) : partials 1 0 y
where x = numerator z
y = denominator z
m = fromInteger $ (42 * abs x) `div` y
x2 = x * x
ty = 2 * y
partials n a b =
let a' = a + n * x2
b' = b + ty
in (a', b') : partials (n + 2) a' b'
arcsin :: FinitePrecision p => Decimal p r -> ExtendedDecimal p
arcsin z = let z' = castUp z
in castDown' $ two * arctan (z' / (one + sqrt (one - z' * z')))
arccos :: FinitePrecision p => Decimal p r -> ExtendedDecimal p
arccos = castDown' . (halfPi -) . arcsin . castUp
seriesPi :: FinitePrecision p => ExtendedDecimal p
seriesPi = castDown' $ 6 * arcsin oneHalf
machinPi :: FinitePrecision p => ExtendedDecimal p
machinPi = castDown' $ 16 * arctan' (1 % 5) - 4 * arctan' (1 % 239)
cfPi :: FinitePrecision p => ExtendedDecimal p
cfPi = pi'
where pi' = continuedFraction m
0 $ (4, 1) : [ (n * n, n * 2 + 1) | n <- [1..] ]
Just p = precision pi'
m = (p `div` 3) * 4
fastPi :: FinitePrecision p => ExtendedDecimal p
fastPi = 3.1415926535897932384626433832795028841971693993751
halfPi :: FinitePrecision p => ExtendedDecimal p
halfPi = castDown' $ pi * oneHalf
quarterPi :: FinitePrecision p => ExtendedDecimal p
quarterPi = castDown' $ pi * oneQuarter
cordic :: FinitePrecision p
=> ExtendedDecimal p -> (ExtendedDecimal p, ExtendedDecimal p)
cordic beta@Num{}
| beta > halfPi = negatePair $ cordic (beta - pi)
| beta < flipSign halfPi = negatePair $ cordic (beta + pi)
| isZero beta = (one, zero)
| otherwise = cordic' beta (one, zero) one angles
where negatePair (x, y) = (flipSign x, flipSign y)
angles = quarterPi : [ arctan' z | let half = 1 % 2
, z <- iterate (* half) half ]
cordic' beta v@(x, y) powerOfTwo ~(angle:angles)
| v' == v = (k * x, k * y)
| otherwise = cordic' beta' v' powerOfTwo' angles
where isNegBeta = isNegative beta
beta' | isNegBeta = beta + angle
| otherwise = beta - angle
factor | isNegBeta = powerOfTwo { sign = Neg }
| otherwise = powerOfTwo
v' = (x - factor * y, factor * x + y)
powerOfTwo' = powerOfTwo * oneHalf
k | p <= 50 = fastK
| otherwise = seriesK
where Just p = precision k
fastK = 0.60725293500888125616944675250492826311239085215009
seriesK = infiniteSeries (*)
[ recip $ sqrt (one + x) | x <- iterate (* oneQuarter) one ]
cordic _ = (qNaN, qNaN)
castUp :: Precision p => Decimal p r -> ExtendedDecimal (PPlus1 (PPlus1 p))
castUp = coerce
castDown' :: Precision p
=> ExtendedDecimal (PPlus1 (PPlus1 p)) -> ExtendedDecimal p
castDown' = cast
castDown :: (Precision p, Rounding r)
=> ExtendedDecimal (PPlus1 (PPlus1 p)) -> Decimal p r
castDown = castRounding . castDown'
instance (FinitePrecision p, Rounding r) => Floating (Decimal p r) where
pi = castRounding pi'
where pi' | p <= 50 = fastPi
| otherwise = cfPi
Just p = precision pi'
exp = castRounding . evalOp . Op.exp
log = castRounding . evalOp . Op.ln
logBase b@Num{} x
| b == ten = castRounding $ evalOp (Op.log10 x)
| x == one = case b `compare` one of
LT -> zero { sign = Neg }
EQ -> qNaN
GT -> zero
| x == b && not (isZero b) = one
logBase b x = evalOp (join $ Op.divide <$> Op.ln x <*> Op.ln b)
x ** y = evalOp (x `Op.power` y)
sqrt = castRounding . evalOp . Op.squareRoot
sin = castDown . snd . cordic . castUp
cos = castDown . fst . cordic . castUp
tan = castDown . uncurry (flip (/)) . cordic . castUp
asin = castRounding . arcsin
acos = castRounding . arccos
atan = castRounding . arctan
sinh x = castDown . evalOp' $
Op.exp x >>= \ex -> two `Op.multiply` ex >>= \tex ->
ex `Op.multiply` ex >>= (`Op.subtract` one) >>= (`Op.divide` tex)
cosh x = castDown . evalOp' $
Op.exp x >>= \ex -> two `Op.multiply` ex >>= \tex ->
ex `Op.multiply` ex >>= (`Op.add` one) >>= (`Op.divide` tex)
tanh x = castDown . evalOp' $
two `Op.multiply` x >>= Op.exp >>= \e2x ->
e2x `Op.subtract` one >>= \e2xm1 -> e2x `Op.add` one >>= (e2xm1 `Op.divide`)
asinh x = castDown . evalOp' $ x `Op.multiply` x >>=
(`Op.add` one) >>= Op.squareRoot >>= (x `Op.add`) >>= Op.ln
acosh x = castDown . evalOp' $ x `Op.multiply` x >>=
(`Op.subtract` one) >>= Op.squareRoot >>= (x `Op.add`) >>= Op.ln
atanh x = castDown . evalOp' $ one `Op.add` x >>= \xp1 ->
one `Op.subtract` x >>= (xp1 `Op.divide`) >>= Op.ln >>= Op.multiply oneHalf
instance (FinitePrecision p, Rounding r) => RealFloat (Decimal p r) where
floatRadix _ = 10
floatDigits x = let Just p = precision x in p
floatRange x = let Just emin = eMin x
Just emax = eMax x
in (fromIntegral emin, fromIntegral emax)
decodeFloat x = case x of
Num { sign = s, coefficient = c, exponent = e }
| c == 0 -> (0, 0)
| otherwise -> (m, n)
where m = signFunc s (fromIntegral $ c * 10^d)
n = fromIntegral e - d
d = floatDigits x - numDigits c
Inf { sign = s } -> (special s 0, maxBound)
NaN { sign = s, signaling = sig, payload = p } -> (special s p, n)
where n = minBound + fromEnum sig
where special :: Sign -> Coefficient -> Integer
special s v = signFunc s (pp + fromIntegral v)
pp = 10 ^ floatDigits x :: Integer
encodeFloat m n = x
where x | am >= pp = special
| otherwise = cast Num { sign = signMatch m
, coefficient = fromInteger am
, exponent = fromIntegral n
}
special
| n == maxBound = Inf { sign = signMatch m }
| n == minBound = qNaN { sign = signMatch m, payload = p }
| otherwise = sNaN { sign = signMatch m, payload = p }
where p = fromInteger (am - pp)
am = abs m :: Integer
pp = 10 ^ floatDigits x :: Integer
isNaN x = case x of
NaN{} -> True
_ -> False
isInfinite x = case x of
Inf{} -> True
_ -> False
isDenormalized = isSubnormal
isNegativeZero x = case x of
Num { sign = Neg, coefficient = 0 } -> True
_ -> False
isIEEE _ = True
instance FinitePrecision p => Bits (Decimal p r) where
x .&. y = evalOp (x `Op.and` y)
x .|. y = evalOp (x `Op.or` y)
x `xor` y = evalOp (x `Op.xor` y)
complement = evalOp . Op.invert
shift x i = evalOp $ Op.shift x (fromIntegral i :: GeneralDecimal)
rotate x i = evalOp $ Op.rotate x (fromIntegral i :: GeneralDecimal)
zeroBits = zero
bit i | i >= 0 && i < finiteBitSize x = x
| otherwise = zeroBits
where x = coerce zero { coefficient = 10 ^ i }
testBit x@Num { sign = Pos, coefficient = c, exponent = 0 } i
| i >= 0 && i < finiteBitSize x = (c `quot` 10 ^ i) `rem` 10 == 1
testBit _ _ = False
bitSizeMaybe = precision
bitSize = finiteBitSize
isSigned _ = False
popCount Num { sign = Pos, coefficient = c, exponent = 0 } = popCount' c 0
where popCount' :: Coefficient -> Int -> Int
popCount' 0 c = c
popCount' x c = case d of
0 -> popCount' x' c
1 -> popCount' x' $! succ c
_ -> 0
where (x', d) = x `quotRem` 10
popCount _ = 0
instance FinitePrecision p => FiniteBits (Decimal p r) where
finiteBitSize x = let Just p = precision x in p
instance NFData (Decimal p r) where
rnf Num { sign = s, coefficient = c, exponent = e } =
rnf s `seq` rnf c `seq` rnf e
rnf Inf { sign = s } = rnf s
rnf NaN { sign = s, signaling = sig, payload = p } =
rnf s `seq` rnf sig `seq` rnf p
zero :: Decimal p r
zero = Num { sign = Pos
, coefficient = 0
, exponent = 0
}
oneHalf :: Decimal p r
oneHalf = zero { coefficient = 5, exponent = -1 }
oneQuarter :: Decimal p r
oneQuarter = zero { coefficient = 25, exponent = -2 }
one :: Decimal p r
one = zero { coefficient = 1 }
two :: Decimal p r
two = zero { coefficient = 2 }
three :: Decimal p r
three = zero { coefficient = 3 }
ten :: Decimal p r
ten = zero { coefficient = 10 }
negativeOne :: Decimal p r
negativeOne = one { sign = Neg }
infinity :: Decimal p r
infinity = Inf { sign = Pos }
qNaN :: Decimal p r
qNaN = NaN { sign = Pos, signaling = False, payload = 0 }
sNaN :: Decimal p r
sNaN = qNaN { signaling = True }
flipSign :: Decimal p r -> Decimal p r
flipSign n = n { sign = negateSign (sign n) }
cast :: (Precision p, Rounding r) => Decimal a b -> Decimal p r
cast = evalOp . roundDecimal
castRounding :: Decimal p a -> Decimal p r
castRounding = coerce
numDigits :: Coefficient -> Int
numDigits x
| x < 10 = 1
| x < 100 = 2
| x < 1000 = 3
| x < 10000 = 4
| x < 100000 = 5
| x < 1000000 = 6
| x < 10000000 = 7
| x < 100000000 = 8
| x < 1000000000 = 9
| otherwise = 9 + numDigits (x `quot` 1000000000)
maxCoefficient :: Precision p => p -> Maybe Coefficient
maxCoefficient p = (\d -> 10 ^ d - 1) <$> precision p
isPositive :: Decimal p r -> Bool
isPositive n = case sign n of
Pos -> True
Neg -> False
isNegative :: Decimal p r -> Bool
isNegative n = case sign n of
Neg -> True
Pos -> False
isFinite :: Decimal p r -> Bool
isFinite Num{} = True
isFinite _ = False
isZero :: Decimal p r -> Bool
isZero Num { coefficient = 0 } = True
isZero _ = False
isNormal :: Precision p => Decimal p r -> Bool
isNormal n
| isFinite n && not (isZero n) = maybe True (adjustedExponent n >=) (eMin n)
| otherwise = False
isSubnormal :: Precision p => Decimal p r -> Bool
isSubnormal n
| isFinite n && not (isZero n) = maybe False (adjustedExponent n <) (eMin n)
| otherwise = False
fromBool :: Bool -> Decimal p r
fromBool False = zero
fromBool True = one
toBool :: Decimal p r -> Bool
toBool Num { coefficient = c } = c /= 0
toBool NaN{} = False
toBool _ = True
fromOrdering :: Ordering -> Decimal p r
fromOrdering LT = negativeOne
fromOrdering EQ = zero
fromOrdering GT = one
eLimit :: Precision p => p -> Maybe Exponent
eLimit = eMax
eTiny :: Precision p => p -> Maybe Exponent
eTiny n = (-) <$> eMin n <*> (fromIntegral . subtract 1 <$> precision n)
eRange :: Precision p => Decimal p r -> Maybe (Exponent, Exponent)
eRange n@Num { coefficient = c } = range <$> eLimit n
where range :: Exponent -> (Exponent, Exponent)
range lim = (-lim - clm1 + 1, lim - clm1)
clength = numDigits c :: Int
clm1 = fromIntegral (clength - 1) :: Exponent
eRange _ = Nothing
adjustedExponent :: Decimal p r -> Exponent
adjustedExponent Num { coefficient = c, exponent = e } =
e + fromIntegral (clength - 1)
where clength = numDigits c :: Int
adjustedExponent _ = error "adjustedExponent: not a finite number"
integralValue :: Decimal a b -> Maybe Integer
integralValue Num { sign = s, coefficient = c, exponent = e }
| e >= 0 = Just (signFunc s $ fromIntegral c * 10^e)
| r == 0 = Just (signFunc s $ fromIntegral q)
where (q, r) = c `quotRem` (10^(-e))
integralValue _ = Nothing