diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2016, Robert Leslie +Copyright (c) 2016-2017, Robert Leslie All rights reserved. diff --git a/TODO b/TODO --- a/TODO +++ b/TODO @@ -1,5 +1,4 @@ -*- Outline -*- * To Do -** instance Floating (Decimal p r) ** instance PrintfArg (Decimal p r) diff --git a/decimal-arithmetic.cabal b/decimal-arithmetic.cabal --- a/decimal-arithmetic.cabal +++ b/decimal-arithmetic.cabal @@ -1,6 +1,6 @@ name: decimal-arithmetic -version: 0.4.0.0 +version: 0.5.0.0 synopsis: An implementation of the General Decimal Arithmetic Specification @@ -16,7 +16,7 @@ license: BSD3 license-file: LICENSE -copyright: © 2016 Robert Leslie +copyright: © 2016–2017 Robert Leslie author: Rob Leslie maintainer: Rob Leslie @@ -38,20 +38,41 @@ hs-source-dirs: src exposed-modules: Numeric.Decimal - Numeric.Decimal.Conversion Numeric.Decimal.Arithmetic + Numeric.Decimal.Conversion + Numeric.Decimal.Encoding Numeric.Decimal.Operation - other-modules: Numeric.Decimal.Number + other-modules: Numeric.Decimal.Exception + Numeric.Decimal.Number Numeric.Decimal.Precision Numeric.Decimal.Rounding - build-depends: base >= 4.7 && < 5 - , mtl + build-depends: base >= 4.8 && < 5 + , binary >= 0.8 && < 0.9 + , binary-bits >= 0.5 && < 0.6 + , deepseq >= 1.4 && < 1.5 + , mtl >= 2.2 && < 2.3 default-language: Haskell2010 default-extensions: Trustworthy other-extensions: FlexibleInstances MultiParamTypeClasses RoleAnnotations + +test-suite spec + type: exitcode-stdio-1.0 + hs-source-dirs: test + main-is: Spec.hs + other-modules: Arbitrary + Numeric.Decimal.EncodingSpec + Numeric.Decimal.NumberSpec + Numeric.Decimal.OperationSpec + build-depends: base + , binary + , decimal-arithmetic + , hspec + , QuickCheck + ghc-options: -threaded -rtsopts -with-rtsopts=-N + default-language: Haskell2010 test-suite doctests type: exitcode-stdio-1.0 diff --git a/src/Numeric/Decimal.hs b/src/Numeric/Decimal.hs --- a/src/Numeric/Decimal.hs +++ b/src/Numeric/Decimal.hs @@ -1,7 +1,7 @@ {-| Module : Numeric.Decimal Description : General arbitrary-precision decimal floating-point number type -Copyright : © 2016 Robert Leslie +Copyright : © 2016–2017 Robert Leslie License : BSD3 Maintainer : rob@mars.org Stability : experimental @@ -38,20 +38,49 @@ representation via 'Show' and 'Read' instances. Note that there may be multiple representations of values that are numerically equal (e.g. 1 and 1.00) which are preserved by this conversion. + +Some decimal numbers also support encoding and decoding specific IEEE 754 +interchange formats via a 'Data.Binary.Binary' instance. -} module Numeric.Decimal ( -- * Usage -- $usage + -- ** Advanced usage + -- $advanced-usage + -- * Arbitrary-precision decimal numbers Decimal , BasicDecimal , ExtendedDecimal , GeneralDecimal + -- ** Number types with defined encodings + -- $encodings + , Decimal32 + , Decimal64 + , Decimal128 + -- ** Precision types - , module Numeric.Decimal.Precision + , Precision + , FinitePrecision + , P1 , P2 , P3 , P4 , P5 , P6 , P7 , P8 , P9 , P10 + , P11, P12, P13, P14, P15, P16, P17, P18, P19, P20 + , P21, P22, P23, P24, P25, P26, P27, P28, P29, P30 + , P31, P32, P33, P34, P35, P36, P37, P38, P39, P40 + , P41, P42, P43, P44, P45, P46, P47, P48, P49, P50 + + , P75, P100, P150, P200, P250, P300, P400, P500, P1000, P2000 + + , PPlus1, PTimes2 + + , PInfinite + + , Pdecimal32 + , Pdecimal64 + , Pdecimal128 + -- ** Rounding types , Rounding @@ -67,23 +96,14 @@ -- * Functions , cast , fromBool + , fromOrdering ) where +import Numeric.Decimal.Encoding import Numeric.Decimal.Number import Numeric.Decimal.Precision import Numeric.Decimal.Rounding --- | A decimal floating point number with 9 digits of precision, rounding half --- up -type BasicDecimal = Decimal P9 RoundHalfUp - --- | A decimal floating point number with selectable precision, rounding half --- even -type ExtendedDecimal p = Decimal p RoundHalfEven - --- | A decimal floating point number with infinite precision -type GeneralDecimal = ExtendedDecimal PInfinite - {- $usage You should choose a decimal number type with appropriate precision and @@ -93,35 +113,52 @@ rounds half up. * 'ExtendedDecimal' is a number type constructor with selectable precision -that rounds half even. For example, @'ExtendedDecimal' 'P34'@ is a number type -with 34 decimal digits of precision. There is a range of ready-made precisions -available, including 'P1' through 'P50' on up to 'P2000' (the IEEE 754 -smallest and basic formats correspond to precisions 'P7', 'P16', or 'P34'). -Alternatively, an arbitrary precision can be constructed through type -application of 'PPlus1' and/or 'PTimes2' to any existing precision. +that rounds half even. For example, @'ExtendedDecimal' 'P15'@ is a number type +with 15 decimal digits of precision. There is a range of ready-made precisions +available, including 'P1' through 'P50' on up to 'P2000'. * 'GeneralDecimal' is a number type with infinite precision. Note that not all operations support numbers with infinite precision. +* 'Decimal32', 'Decimal64', and 'Decimal128' are specialized number types with +'Data.Binary.Binary' instances that implement the /decimal32/, /decimal64/, +and /decimal128/ interchange format encodings described in IEEE 754-2008. +These types have precisions of 7, 16, and 34 decimal digits, respectively, and +round half even. + * The most versatile 'Decimal' type constructor is parameterized by both a precision and a rounding algorithm. For example, @'Decimal' 'P20' 'RoundDown'@ is a number type with 20 decimal digits of precision that rounds down (truncates). Several 'Rounding' algorithms are available to choose from. -It is suggested to create an alias for the type of numbers you wish to support -in your application. For example: - -> type Number = ExtendedDecimal P16 - -A decimal number type may be used in a @default@ declaration, possibly -replacing 'Double' and/or 'Integer'. For example: +A decimal number type may be used in a @default@ declaration, for example +replacing 'Double': > default (Integer, BasicDecimal) +-} -== Advanced usage +{- $advanced-usage -Additional operations and control beyond what is provided by the basic numeric -type classes are available through the use of "Numeric.Decimal.Arithmetic" and +Additional operations and control beyond what is provided by the standard type +classes are available through the use of "Numeric.Decimal.Arithmetic" and "Numeric.Decimal.Operation". Advanced string conversion is also available through "Numeric.Decimal.Conversion". + +Arbitrary precisions can be constructed through type application of 'PPlus1' +and/or 'PTimes2' to any existing precision. + +It is possible to create arbitrary width interchange format encodings with the +help of "Numeric.Decimal.Encoding". +-} + +{- $encodings + +These decimal number types have a 'Data.Binary.Binary' instance that +implements a specific interchange format encoding described in IEEE +754-2008. See "Numeric.Decimal.Encoding" for further details, including the +ability to create additional formats of arbitrary width. + +Alternative rounding algorithms can be used through the more general 'Decimal' +type constructor and the special precision types 'Pdecimal32', 'Pdecimal64', +or 'Pdecimal128', e.g. @'Decimal' 'Pdecimal64' 'RoundCeiling'@. -} diff --git a/src/Numeric/Decimal/Arithmetic.hs b/src/Numeric/Decimal/Arithmetic.hs --- a/src/Numeric/Decimal/Arithmetic.hs +++ b/src/Numeric/Decimal/Arithmetic.hs @@ -71,8 +71,8 @@ -- arithmetic computation and manipulate its 'Context'. -- | A context for decimal arithmetic, carrying signal flags, trap enabler --- state, and a trap handler, parameterized by precision @p@ and rounding --- algorithm @r@ +-- state, and a trap handler, parameterized by precision @p@ and rounding mode +-- @r@ data Context p r = Context { flags :: Signals -- ^ The current signal flags of the context @@ -105,9 +105,8 @@ disabled = [Inexact, Rounded, Subnormal] -- | Return a new context with all signal flags cleared, all traps disabled --- (IEEE 854 §7), using selectable precision (the IEEE 754 smallest and basic --- formats correspond to precisions 'P7', 'P16', or 'P34'), and rounding half --- even (IEEE 754 §4.3.3). +-- (IEEE 854 §7), using selectable precision, and rounding half even (IEEE 754 +-- §4.3.3). extendedDefaultContext :: Context p RoundHalfEven extendedDefaultContext = newContext @@ -123,7 +122,7 @@ deriving Show -- | A decimal arithmetic monad parameterized by the precision @p@ and --- rounding algorithm @r@ +-- rounding mode @r@ newtype Arith p r a = Arith (ExceptT (Exception p r) (State (Context p r)) a) @@ -164,8 +163,8 @@ evalArith (Arith e) = evalState (runExceptT e) -- | Perform a subcomputation using a different precision and/or rounding --- algorithm. The subcomputation is evaluated within a new context with all --- flags cleared and all traps disabled. Any flags set in the context of the +-- mode. The subcomputation is evaluated within a new context with all flags +-- cleared and all traps disabled. Any flags set in the context of the -- subcomputation are ignored, but if an exception is returned it will be -- re-raised within the current context. subArith :: Arith a b (Decimal a b) -> Arith p r (Decimal a b) @@ -181,7 +180,7 @@ where getPrecision' :: Precision p => p -> Arith p r (Maybe Int) getPrecision' = return . precision --- | Return the rounding algorithm of the arithmetic context. +-- | Return the rounding mode of the arithmetic context. getRounding :: Rounding r => Arith p r RoundingAlgorithm getRounding = getRounding' undefined where getRounding' :: Rounding r => r -> Arith p r RoundingAlgorithm diff --git a/src/Numeric/Decimal/Conversion.hs b/src/Numeric/Decimal/Conversion.hs --- a/src/Numeric/Decimal/Conversion.hs +++ b/src/Numeric/Decimal/Conversion.hs @@ -244,15 +244,18 @@ showString (replicate (fromIntegral $ -e - cl) '0') . showString cs - Inf { } -> showString "Infinity" - QNaN { payload = p } -> showString "NaN" . diag p - SNaN { payload = p } -> showString "sNaN" . diag p + Inf { } -> showString "Infinity" + NaN { signaling = s, payload = p } -> sig s . showString "NaN" . diag p where signStr :: ShowS signStr = case sign num of Pos -> id Neg -> showChar '-' + sig :: Bool -> ShowS + sig False = id + sig True = showChar 's' + diag :: Payload -> ShowS diag 0 = id diag d = shows d @@ -295,34 +298,26 @@ , exponent = e - fromIntegral (length fracDigits) } - digitsWithOptionalPoint = fractionalDigits <|> wholeDigits - - fractionalDigits = do - char '.' - fracDigits <- many1 parseDigit - return $ \e -> - Num { sign = Pos - , coefficient = readDigits fracDigits - , exponent = e - fromIntegral (length fracDigits) - } - - wholeDigits = do + digitsWithOptionalPoint = do + fractional <- option False (char '.' *> pure True) digits <- many1 parseDigit + let offset | fractional = fromIntegral (length digits) + | otherwise = 0 return $ \e -> Num { sign = Pos , coefficient = readDigits digits - , exponent = e + , exponent = e - offset } parseExponentPart :: ReadP Exponent parseExponentPart = do parseString "E" - parseSign negate <*> (readDigits <$> many1 parseDigit) + parseSign negate <*> (fromIntegral . readDigits <$> many1 parseDigit) parseInfinity :: ReadP (Decimal p r) parseInfinity = do parseString "Inf" optional $ parseString "inity" - return Inf { sign = Pos } + return infinity parseNaN :: ReadP (Decimal p r) parseNaN = parseQNaN <|> parseSNaN @@ -330,13 +325,13 @@ parseQNaN :: ReadP (Decimal p r) parseQNaN = do p <- parseNaNPayload - return QNaN { sign = Pos, payload = p } + return qNaN { payload = p } parseSNaN :: ReadP (Decimal p r) parseSNaN = do parseString "s" p <- parseNaNPayload - return SNaN { sign = Pos, payload = p } + return sNaN { payload = p } parseNaNPayload :: ReadP Payload parseNaNPayload = do @@ -349,7 +344,7 @@ parseString :: String -> ReadP () parseString = mapM_ $ \c -> char (toLower c) <|> char (toUpper c) - readDigits :: Num c => [Int] -> c + readDigits :: [Int] -> Coefficient readDigits = foldl' (\a b -> a * 10 + fromIntegral b) 0 {- $doctest-toNumber diff --git a/src/Numeric/Decimal/Encoding.hs b/src/Numeric/Decimal/Encoding.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Decimal/Encoding.hs @@ -0,0 +1,330 @@ + +{-# LANGUAGE FlexibleInstances #-} + +-- | This module implements the decimal interchange format encodings described +-- in IEEE 754-2008, including the /decimal32/, /decimal64/, and /decimal128/ +-- formats, as well as arbitrary width /decimal{k}/ formats through the use of +-- 'Format' with 'KPlus32' and\/or 'KTimes2'. For example, to use a +-- /decimal96/ format: +-- +-- > type Decimal96 = ExtendedDecimal (Format (KPlus32 K64) DecimalCoefficient) +-- +-- Currently only a decimal encoding of coefficients is implemented, but a +-- binary encoding may be added in the future. +module Numeric.Decimal.Encoding ( + -- * Primary convenience types + Decimal32 + , Decimal64 + , Decimal128 + + -- ** Precision types + , Pdecimal32 + , Pdecimal64 + , Pdecimal128 + + -- * Interchange format types + , Format + + -- ** Format parameters + , Parameters + , K32 + , K64 + , K128 + , KPlus32 + , KTimes2 + + -- ** Coefficient encodings + -- , CoefficientEncoding + , DecimalCoefficient + , BinaryCoefficient + ) where + +import Prelude hiding (exponent) + +import Data.Binary (Binary(get, put), Get) +import Data.Binary.Bits.Get (BitGet, getBool, getWord8, getWord16be, runBitGet) +import Data.Binary.Bits.Put (BitPut, putBool, putWord8, putWord16be, runBitPut) +import Data.Bits (bit, shiftL, shiftR, testBit, (.&.), (.|.)) +import Data.Word (Word8, Word16) + +import Numeric.Decimal.Number +import Numeric.Decimal.Precision + +-- Decimal number types + +-- | A decimal floating point number with 7 digits of precision, rounding half +-- even, and a 32-bit encoded representation using the /decimal32/ interchange +-- format (with a decimal encoding for the coefficient) +type Decimal32 = ExtendedDecimal Pdecimal32 + +-- | A decimal floating point number with 16 digits of precision, rounding +-- half even, and a 64-bit encoded representation using the /decimal64/ +-- interchange format (with a decimal encoding for the coefficient) +type Decimal64 = ExtendedDecimal Pdecimal64 + +-- | A decimal floating point number with 34 digits of precision, rounding +-- half even, and a 128-bit encoded representation using the /decimal128/ +-- interchange format (with a decimal encoding for the coefficient) +type Decimal128 = ExtendedDecimal Pdecimal128 + +-- Precision types + +-- | A type with 'Precision' instance specifying /decimal32/ interchange +-- format parameters (using a decimal encoding for the coefficient) having an +-- effective precision of 7 decimal digits +type Pdecimal32 = Format K32 DecimalCoefficient + +-- | A type with 'Precision' instance specifying /decimal64/ interchange +-- format parameters (using a decimal encoding for the coefficient) having an +-- effective precision of 16 decimal digits +type Pdecimal64 = Format K64 DecimalCoefficient + +-- | A type with 'Precision' instance specifying /decimal128/ interchange +-- format parameters (using a decimal encoding for the coefficient) having an +-- effective precision of 34 decimal digits +type Pdecimal128 = Format K128 DecimalCoefficient + +-- Format parameters + +-- | Interchange format parameters used to define an encoding and derive the +-- format's /precision/ and E/max/ +class Parameters k where + -- | /k//32, the primary format parameter defining the encoding width as a + -- multiple of 32 bits + paramK32 :: k -> Int + +-- | /p/, precision in digits +paramP :: Parameters k => k -> Int +paramP k = 9 * paramK32 k - 2 + +-- | /emax/ +paramEmax :: Parameters k => k -> Exponent +paramEmax k = 3 * 2^(paramK32 k * 2 + 3) + +-- | /bias/, /E/ − /q/ +paramBias :: Parameters k => k -> Exponent +paramBias k = paramEmax k + fromIntegral (paramP k - 2) + +-- | /w/, combination field width in bits − 5 +paramW :: Parameters k => k -> Int +paramW k = paramK32 k * 2 + 4 + +-- | /t//10, trailing significand field width in 10-bit multiples +paramT10 :: Parameters k => k -> Int +paramT10 k = 3 * paramK32 k - 1 + +-- | Parameters for the /decimal32/ interchange format +data K32 +instance Parameters K32 where + paramK32 _ = 1 + +-- | Parameters for the /decimal64/ interchange format +type K64 = KPlus32 K32 + +-- | Parameters for the /decimal128/ interchange format +type K128 = KTimes2 K64 + +-- | Parameters for a /decimal{@k@ + 32}/ interchange format +data KPlus32 k +instance Parameters k => Parameters (KPlus32 k) where + paramK32 t = paramK32 (minus32 t) + 1 + where minus32 :: KPlus32 k -> k + minus32 = undefined + +-- | Parameters for a /decimal{@k@ × 2}/ interchange format +data KTimes2 k +instance Parameters k => Parameters (KTimes2 k) where + paramK32 t = paramK32 (div2 t) * 2 + where div2 :: KTimes2 k -> k + div2 = undefined + +-- | A class encapsulating coefficient encodings +class CoefficientEncoding c + +-- | Specify a decimal encoding for the coefficient. +data DecimalCoefficient +instance CoefficientEncoding DecimalCoefficient + +-- | Specify a binary encoding for the coefficient (currently unimplemented). +data BinaryCoefficient +instance CoefficientEncoding BinaryCoefficient + +-- | A type (with a 'Precision' instance) for specifying interchange format +-- parameters @k@ and coefficient encoding @c@ +data Format k c + +formatK :: Format k c -> k +formatK = undefined + +-- | This 'Precision' instance automatically computes the /precision/ and +-- E/max/ of decimal numbers that use this format. +instance Parameters k => Precision (Format k c) where + precision = Just . paramP . formatK + eMax = Just . paramEmax . formatK + +instance Parameters k => FinitePrecision (Format k c) + +-- | A 'Binary' instance is defined for interchange formats for which a +-- 'Parameters' instance exists, and covers particularly the 'Decimal32', +-- 'Decimal64', and 'Decimal128' types. +instance Parameters k => Binary (Decimal (Format k DecimalCoefficient) r) where + put d = runBitPut $ putDecimal (paramW k) (paramT10 k) (paramBias k) d + where k = formatK (decimalFormat d) + + decimalFormat :: Decimal (Format k c) r -> Format k c + decimalFormat = undefined + + get = result + where result = runBitGet $ getDecimal (paramW k) (paramT10 k) (paramBias k) + k = formatK (getDecimalFormat result) + + getDecimalFormat :: Get (Decimal (Format k c) r) -> Format k c + getDecimalFormat = undefined + +-- Densely Packed Decimal + +dpd2bcd :: Word16 -> (Word8, Word8, Word8) +dpd2bcd dpd = case mask 0 0xe of + 0xe -> case mask 4 0x6 of + 0x6 -> ( c7, f4, i0) + 0x4 -> (a9b8c7, f4, i0) + 0x2 -> ( c7, d9e8f4, i0) + _ -> ( c7, f4, g9h8i0) + 0xc -> ( c7, d6e5f4, g9h8i0) + 0xa -> (a9b8c7, f4, g6h5i0) + 0x8 -> (a9b8c7, d6e5f4, i0) + _ -> (a9b8c7, d6e5f4, g2h1i0) + + where a9b8c7 = mask 7 7 + d6e5f4 = mask 4 7 + d9e8f4 = mask 7 6 .|. mask 4 1 + g2h1i0 = mask 0 7 + g6h5i0 = mask 4 6 .|. mask 0 1 + g9h8i0 = mask 7 6 .|. mask 0 1 + i0 = 8 .|. mask 0 1 + f4 = 8 .|. mask 4 1 + c7 = 8 .|. mask 7 1 + + mask :: Int -> Word8 -> Word8 + mask s m = fromIntegral (shiftR dpd s) .&. m + +bcd2dpd :: Word8 -> Word8 -> Word8 -> Word16 +bcd2dpd d2 d1 d0 = case (d2 < 8, d1 < 8, d0 < 8) of + (True , True , True ) -> a9b8c7 .|. d6e5f4 .|. g2h1i0 + (True , True , False) -> a9b8c7 .|. d6e5f4 .|. 0x08 .|. i0 + (True , False, True ) -> a9b8c7 .|. g6h5 .|. f4 .|. 0x0a .|. i0 + (False, True , True ) -> g9h8 .|. c7 .|. d6e5f4 .|. 0x0c .|. i0 + (False, False, True ) -> g9h8 .|. c7 .|. f4 .|. 0x0e .|. i0 + (False, True , False) -> d9e8 .|. c7 .|. f4 .|. 0x2e .|. i0 + (True , False, False) -> a9b8c7 .|. f4 .|. 0x4e .|. i0 + (False, False, False) -> c7 .|. f4 .|. 0x6e .|. i0 + + where a9b8c7 = isolate d2 7 7 + c7 = isolate d2 1 7 + d6e5f4 = isolate d1 7 4 + d9e8 = isolate d1 6 7 + f4 = isolate d1 1 4 + g2h1i0 = isolate d0 7 0 + g6h5 = isolate d0 6 4 + g9h8 = isolate d0 6 7 + i0 = isolate d0 1 0 + + isolate :: Word8 -> Word8 -> Int -> Word16 + isolate d m = shiftL (fromIntegral $ d .&. m) + +-- Low-level encoding/decoding + +data CombinationField = Finite { exponentMSBs :: Word8 + , coefficientMSD :: Word8 } + | Infinity + | NotANumber + +getCommon :: BitGet (Sign, CombinationField) +getCommon = do + sign <- toEnum . fromEnum <$> getBool + bits <- getWord8 5 + let cf = case bits of + 0x1e -> Infinity + 0x1f -> NotANumber + _ -> let ab = shiftR bits 3 in case ab of + 0x03 -> Finite { exponentMSBs = shiftR bits 1 .&. 0x03 + , coefficientMSD = 0x08 .|. (bits .&. 0x01) + } + _ -> Finite { exponentMSBs = ab + , coefficientMSD = bits .&. 0x07 + } + return (sign, cf) + +putCommon :: Sign -> CombinationField -> BitPut () +putCommon sign cf = do + putBool (toEnum . fromEnum $ sign) + let bits = case cf of + Finite { exponentMSBs = msbs, coefficientMSD = msd } + | msd < 8 -> shiftL msbs 3 .|. msd + | otherwise -> 0x18 .|. shiftL msbs 1 .|. (msd .&. 0x01) + Infinity -> 0x1e + NotANumber -> 0x1f + putWord8 5 bits + +getCoefficient :: CombinationField -> Int -> BitGet Coefficient +getCoefficient = getCoefficient' . getMSD + + where getCoefficient' :: Coefficient -> Int -> BitGet Coefficient + getCoefficient' ic 0 = return ic + getCoefficient' ic n = do + (a, b, c) <- dpd2bcd <$> getWord16be 10 + let v = fromIntegral a * 100 + fromIntegral b * 10 + fromIntegral c + getCoefficient' (ic * 1000 + v) (pred n) + + getMSD :: CombinationField -> Coefficient + getMSD Finite { coefficientMSD = msd } = fromIntegral msd + getMSD _ = 0 + +getDecimal :: Int -> Int -> Exponent -> BitGet (Decimal p r) +getDecimal ecbits cclen bias = do + (sign, cf) <- getCommon + ec <- getWord16be ecbits + coefficient <- getCoefficient cf cclen + return $ case cf of + Finite { exponentMSBs = msbs } -> + let ee = shiftL (fromIntegral msbs) ecbits .|. fromIntegral ec + in Num { sign = sign, coefficient = coefficient, exponent = ee - bias } + Infinity -> Inf { sign = sign } + NotANumber -> + let s = testBit ec (ecbits - 1) + in NaN { sign = sign, signaling = s, payload = coefficient } + +putDecimal :: Int -> Int -> Exponent -> Decimal p r -> BitPut () +putDecimal ecbits cclen bias x = do + let msd : cc = digits x + (cf, ee) = case x of + Num { exponent = e } -> + let cf = Finite { exponentMSBs = fromIntegral (shiftR ee ecbits) + , coefficientMSD = msd + } + in (cf, fromIntegral $ e + bias) + Inf{} -> (Infinity, 0) + NaN { signaling = s } -> (NotANumber, if s then bit (ecbits - 1) else 0) + putCommon (sign x) cf + putWord16be ecbits (ee .&. (bit ecbits - 1)) + putDigits cc + + where digits :: Decimal p r -> [Word8] + digits x = let ds = case x of + Num { coefficient = c } -> digits' c + NaN { payload = p } -> digits' p + Inf { } -> [] + in replicate (1 + cclen * 3 - length ds) 0 ++ ds + + digits' :: Coefficient -> [Word8] + digits' = go [] + where go ds 0 = ds + go ds c = let (q, r) = c `quotRem` 10 + in go (fromIntegral r : ds) q + + putDigits :: [Word8] -> BitPut () + putDigits (a : b : c : rest) = do + putWord16be 10 (bcd2dpd a b c) + putDigits rest + putDigits [] = return () + putDigits _ = error "putDigits: invalid # digits" diff --git a/src/Numeric/Decimal/Exception.hs b/src/Numeric/Decimal/Exception.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Decimal/Exception.hs @@ -0,0 +1,242 @@ + +module Numeric.Decimal.Exception ( + -- * Exceptional conditions + clamped + , conversionSyntax + , divisionByZero + , divisionImpossible + , divisionUndefined + , inexact + , insufficientStorage + , invalidContext + , invalidOperation + , overflow + , rounded + , subnormal + , underflow + ) where + +import Prelude hiding (exponent) + +import Numeric.Decimal.Arithmetic +import Numeric.Decimal.Number +import Numeric.Decimal.Precision +import Numeric.Decimal.Rounding + +-- | This occurs and signals 'Clamped' if the exponent of a result has been +-- altered in order to fit the constraints of a specific concrete +-- representation. This may occur when the exponent of a zero result would be +-- outside the bounds of a representation, or (in the IEEE 754 interchange +-- formats) when a large normal number would have an encoded exponent that +-- cannot be represented. In this latter case, the exponent is reduced to fit +-- and the corresponding number of zero digits are appended to the coefficient +-- (“fold-down”). The condition always occurs when a subnormal value rounds to +-- zero. +clamped :: Decimal p r -> Arith p r (Decimal p r) +clamped = raiseSignal Clamped + +-- | This occurs and signals 'InvalidOperation' if a string is being converted +-- to a number and it does not conform to the numeric string syntax. The +-- result is @[0,qNaN]@. +conversionSyntax :: Arith p r (Decimal p r) +conversionSyntax = raiseSignal InvalidOperation qNaN + +-- | This occurs and signals 'DivisionByZero' if division of a finite number +-- by zero was attempted (during a 'Numeric.Decimal.Operation.divideInteger' +-- or 'Numeric.Decimal.Operation.divide' operation, or a +-- 'Numeric.Decimal.Operation.power' operation with negative right-hand +-- operand), and the dividend was not zero. +-- +-- The result of the operation is @[@/sign/@,inf]@, where /sign/ is the +-- exclusive or of the signs of the operands for divide, or is 1 for an odd +-- power of −0, for power. +divisionByZero :: Decimal p r -> Arith p r (Decimal p r) +divisionByZero = raiseSignal DivisionByZero + +-- | This occurs and signals 'InvalidOperation' if the integer result of a +-- 'Numeric.Decimal.Operation.divideInteger' or +-- 'Numeric.Decimal.Operation.remainder' operation had too many digits (would +-- be longer than /precision/). The result is @[0,qNaN]@. +divisionImpossible :: Arith p r (Decimal p r) +divisionImpossible = raiseSignal InvalidOperation qNaN + +-- | This occurs and signals 'InvalidOperation' if division by zero was +-- attempted (during a 'Numeric.Decimal.Operation.divideInteger', +-- 'Numeric.Decimal.Operation.divide', or +-- 'Numeric.Decimal.Operation.remainder' operation), and the dividend is also +-- zero. The result is @[0,qNaN]@. +divisionUndefined :: Arith p r (Decimal p r) +divisionUndefined = raiseSignal InvalidOperation qNaN + +-- | This occurs and signals 'Inexact' whenever the result of an operation is +-- not exact (that is, it needed to be rounded and any discarded digits were +-- non-zero), or if an overflow or underflow condition occurs. The result in +-- all cases is unchanged. +-- +-- The 'Inexact' signal may be tested (or trapped) to determine if a given +-- operation (or sequence of operations) was inexact. +inexact :: Decimal p r -> Arith p r (Decimal p r) +inexact = raiseSignal Inexact + +-- | For many implementations, storage is needed for calculations and +-- intermediate results, and on occasion an arithmetic operation may fail due +-- to lack of storage. This is considered an operating environment error, +-- which can be either be handled as appropriate for the environment, or +-- treated as an Invalid operation condition. The result is @[0,qNaN]@. +insufficientStorage :: Arith p r (Decimal p r) +insufficientStorage = invalidOperation qNaN + +-- | This occurs and signals 'InvalidOperation' if an invalid context was +-- detected during an operation. This can occur if contexts are not checked on +-- creation and either the /precision/ exceeds the capability of the +-- underlying concrete representation or an unknown or unsupported /rounding/ +-- was specified. These aspects of the context need only be checked when the +-- values are required to be used. The result is @[0,qNaN]@. +invalidContext :: Arith p r (Decimal p r) +invalidContext = raiseSignal InvalidOperation qNaN + +-- | This occurs and signals 'InvalidOperation' if: +-- +-- * an operand to an operation is @[s,sNaN]@ or @[s,sNaN,d]@ (any /signaling/ +-- NaN) +-- +-- * an attempt is made to add @[0,inf]@ to @[1,inf]@ during an addition or +-- subtraction operation +-- +-- * an attempt is made to multiply 0 by @[0,inf]@ or @[1,inf]@ +-- +-- * an attempt is made to divide either @[0,inf]@ or @[1,inf]@ by either +-- @[0,inf]@ or @[1,inf]@ +-- +-- * the divisor for a remainder operation is zero +-- +-- * the dividend for a remainder operation is either @[0,inf]@ or @[1,inf]@ +-- +-- * either operand of the 'Numeric.Decimal.Operation.quantize' operation is +-- infinite, or the result of a 'Numeric.Decimal.Operation.quantize' operation +-- would require greater precision than is available +-- +-- * the operand of the 'Numeric.Decimal.Operation.ln' or the +-- 'Numeric.Decimal.Operation.log10' operation is less than zero +-- +-- * the operand of the 'Numeric.Decimal.Operation.squareRoot' operation has a +-- /sign/ of 1 and a non-zero /coefficient/ +-- +-- * both operands of the 'Numeric.Decimal.Operation.power' operation are +-- zero, or if the left-hand operand is less than zero and the right-hand +-- operand does not have an integral value or is infinite +-- +-- * an operand is invalid; for example, certain values of concrete +-- representations may not correspond to numbers — an implementation is +-- permitted (but is not required) to detect these invalid values and raise +-- this condition. +-- +-- The result of the operation after any of these invalid operations is +-- @[0,qNaN]@ except when the cause is a signaling NaN, in which case the +-- result is @[s,qNaN]@ or @[s,qNaN,d]@ where the sign and diagnostic are +-- copied from the signaling NaN. +invalidOperation :: Decimal a b -> Arith p r (Decimal p r) +invalidOperation n = raiseSignal InvalidOperation $ case n of + NaN { signaling = True } -> n { signaling = False } + _ -> qNaN + +-- | This occurs and signals 'Overflow' if the /adjusted exponent/ of a result +-- (from a conversion or from an operation that is not an attempt to divide by +-- zero), after rounding, would be greater than the largest value that can be +-- handled by the implementation (the value E/max/). +-- +-- The result depends on the rounding mode: +-- +-- * For 'RoundHalfUp' and 'RoundHalfEven' (and for 'RoundHalfDown' and +-- 'RoundUp', if implemented), the result of the operation is [sign,@inf@], +-- where /sign/ is the sign of the intermediate result. +-- +-- * For 'RoundDown', (and 'Round05Up', if implemented), the result is the +-- largest finite number that can be represented in the current /precision/, +-- with the sign of the intermediate result. +-- +-- * For 'RoundCeiling', the result is the same as for 'RoundDown' if the sign +-- of the intermediate result is 1, or is @[0,inf]@ otherwise. +-- +-- * For 'RoundFloor', the result is the same as for 'RoundDown' if the sign +-- of the intermediate result is 0, or is @[1,inf]@ otherwise. +-- +-- In all cases, 'inexact' and 'rounded' will also be raised. +-- +-- Note: IEEE 854 §7.3 requires that the result delivered to a trap handler be +-- different, depending on whether the overflow was the result of a conversion +-- or of an arithmetic operation. This specification deviates from IEEE 854 in +-- this respect; however, an implementation could comply with IEEE 854 by +-- providing a separate mechanism for the special result to a trap +-- handler. IEEE 754 has no such requirement. +overflow :: (Precision p, Rounding r) => Decimal p r -> Arith p r (Decimal p r) +overflow ir = result >>= raiseSignal Overflow >>= inexact >>= rounded + + where result :: (Precision p, Rounding r) => Arith p r (Decimal p r) + result = getRounding >>= \r -> case r of + RoundHalfUp -> signedInfinity + RoundHalfEven -> signedInfinity + RoundHalfDown -> signedInfinity + RoundUp -> signedInfinity + + RoundDown -> largestFinite + Round05Up -> largestFinite + + RoundCeiling -> case sign ir of + Neg -> largestFinite + Pos -> return infinity + + RoundFloor -> case sign ir of + Pos -> largestFinite + Neg -> return infinity { sign = Neg } + + signedInfinity :: Arith p r (Decimal p r) + signedInfinity = return infinity { sign = sign ir } + + largestFinite :: Precision p => Arith p r (Decimal p r) + largestFinite = getPrecision >>= \p -> + let x = Num { sign = sign ir + , coefficient = undefined -- 10^p - 1 + , exponent = undefined -- eMax x + } + in return x + +-- | This occurs and signals 'Rounded' whenever the result of an operation is +-- rounded (that is, some zero or non-zero digits were discarded from the +-- coefficient), or if an overflow or underflow condition occurs. The result +-- in all cases is unchanged. +-- +-- The 'Rounded' signal may be tested (or trapped) to determine if a given +-- operation (or sequence of operations) caused a loss of precision. +rounded :: Decimal p r -> Arith p r (Decimal p r) +rounded = raiseSignal Rounded + +-- | This occurs and signals 'Subnormal' whenever the result of a conversion +-- or operation is subnormal (that is, its adjusted exponent is less than +-- E/min/, before any rounding). The result in all cases is unchanged. +-- +-- The 'Subnormal' signal may be tested (or trapped) to determine if a given +-- or operation (or sequence of operations) yielded a subnormal result. +subnormal :: Decimal p r -> Arith p r (Decimal p r) +subnormal = raiseSignal Subnormal + +-- | This occurs and signals 'Underflow' if a result is inexact and the +-- /adjusted exponent/ of the result would be smaller (more negative) than the +-- smallest value that can be handled by the implementation (the value +-- E/min/). That is, the result is both inexact and subnormal. +-- +-- The result after an underflow will be a subnormal number rounded, if +-- necessary, so that its exponent is not less than E/tiny/. This may result +-- in 0 with the sign of the intermediate result and an exponent of E/tiny/. +-- +-- In all cases, 'inexact', 'rounded', and 'subnormal' will also be raised. +-- +-- Note: IEEE 854 §7.4 requires that the result delivered to a trap handler be +-- different, depending on whether the underflow was the result of a +-- conversion or of an arithmetic operation. This specification deviates from +-- IEEE 854 in this respect; however, an implementation could comply with IEEE +-- 854 by providing a separate mechanism for the result to a trap +-- handler. IEEE 754 has no such requirement. +underflow :: Decimal p r -> Arith p r (Decimal p r) +underflow x = undefined >>= raiseSignal Underflow >>= + subnormal >>= inexact >>= rounded diff --git a/src/Numeric/Decimal/Exception.hs-boot b/src/Numeric/Decimal/Exception.hs-boot new file mode 100644 --- /dev/null +++ b/src/Numeric/Decimal/Exception.hs-boot @@ -0,0 +1,11 @@ + +module Numeric.Decimal.Exception ( + inexact + , rounded + ) where + +import {-# SOURCE #-} Numeric.Decimal.Arithmetic +import {-# SOURCE #-} Numeric.Decimal.Number + +inexact :: Decimal p r -> Arith p r (Decimal p r) +rounded :: Decimal p r -> Arith p r (Decimal p r) diff --git a/src/Numeric/Decimal/Number.hs b/src/Numeric/Decimal/Number.hs --- a/src/Numeric/Decimal/Number.hs +++ b/src/Numeric/Decimal/Number.hs @@ -11,6 +11,10 @@ , Payload , Decimal(..) + , BasicDecimal + , ExtendedDecimal + , GeneralDecimal + , zero , oneHalf , one @@ -23,7 +27,9 @@ , flipSign , cast + , fromBool + , fromOrdering , isPositive , isNegative @@ -38,6 +44,7 @@ import Prelude hiding (exponent) +import Control.DeepSeq (NFData(..)) import Control.Monad (join) import Data.Bits (Bits(..), FiniteBits(..)) import Data.Char (isSpace) @@ -55,14 +62,13 @@ import qualified GHC.Real -{- $setup ->>> :load Harness --} - data Sign = Pos -- ^ Positive or non-negative | Neg -- ^ Negative deriving (Eq, Enum) +instance NFData Sign where + rnf s = s `seq` () + negateSign :: Sign -> Sign negateSign Pos = Neg negateSign Neg = Pos @@ -83,25 +89,34 @@ _ -> Pos type Coefficient = Natural -type Exponent = Int +type Exponent = Integer type Payload = Coefficient -- | A decimal floating point number with selectable precision and rounding -- algorithm data Decimal p r - = Num { sign :: Sign - , coefficient :: Coefficient - , exponent :: Exponent - } - | Inf { sign :: Sign - } - | QNaN { sign :: Sign - , payload :: Payload - } - | SNaN { sign :: Sign - , payload :: Payload - } + = Num { sign :: Sign + , coefficient :: Coefficient + , exponent :: Exponent + } + | Inf { sign :: Sign + } + | NaN { sign :: Sign + , signaling :: Bool + , payload :: Payload + } +-- | A decimal floating point number with 9 digits of precision, rounding half +-- up +type BasicDecimal = Decimal P9 RoundHalfUp + +-- | A decimal floating point number with selectable precision, rounding half +-- even +type ExtendedDecimal p = Decimal p RoundHalfEven + +-- | A decimal floating point number with infinite precision +type GeneralDecimal = ExtendedDecimal PInfinite + -- | The 'Show' instance uses the 'toScientificString' operation from -- "Numeric.Decimal.Conversion". instance Show (Decimal p r) where @@ -115,67 +130,61 @@ | (n, s) <- readParen False (readP_to_S toNumber . dropWhile isSpace) str ] -{- $doctest-Read ->>> fmap toRep (read "Just 123" :: Maybe GeneralDecimal) -Just (N (0,123,0)) - ->>> fmap toRep (read "Just (-12.0)" :: Maybe GeneralDecimal) -Just (N (1,120,-1)) --} +decimalPrecision :: Decimal p r -> p +decimalPrecision = undefined instance Precision p => Precision (Decimal p r) where precision = precision . decimalPrecision - where decimalPrecision :: Decimal p r -> p - decimalPrecision = undefined + eMax = eMax . decimalPrecision + eMin = eMin . decimalPrecision -evalOp :: Arith p r (Decimal p r) -> Decimal p r -evalOp op = either exceptionResult id $ evalArith op newContext +-- This assumes the arithmetic operation does not trap any signals, which +-- could result in an exception being thrown (and returned in a Left value). +evalOp :: Arith p r a -> a +evalOp op = let Right r = evalArith op newContext in r -type GeneralDecimal = Decimal PInfinite RoundHalfEven +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) + +-- | Note that NaN values are not equal to any value, including other NaNs. instance Eq (Decimal p r) where - x == y = case evalOp (x `Op.compare` y) :: GeneralDecimal of - Num { coefficient = 0 } -> True - _ -> False + x == y = case compareDecimal x y of + Right EQ -> True + _ -> False -instance (Precision p, Rounding r) => Ord (Decimal p r) where - x `compare` y = case evalOp (x `Op.compare` y) :: GeneralDecimal of - Num { coefficient = 0 } -> EQ - Num { sign = Neg } -> LT - Num { sign = Pos } -> GT - _ -> GT -- match Prelude behavior for NaN +-- | Unlike the instances for 'Float' and 'Double', the 'compare' method in +-- this instance uses a total ordering over all possible values. Note that +-- @'compare' x y == 'EQ'@ does not imply @x == y@ (and similarly for 'LT' and +-- 'GT') in the cases where @x@ or @y@ are NaN values. +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 evalOp (x `Op.compare` y) :: GeneralDecimal of - Num { sign = Neg } -> True - _ -> False + x < y = case compareDecimal x y of + Right LT -> True + _ -> False - x <= y = case evalOp (x `Op.compare` y) :: GeneralDecimal of - Num { sign = Neg } -> True - Num { coefficient = 0 } -> True - _ -> False + x <= y = case compareDecimal x y of + Right LT -> True + Right EQ -> True + _ -> False - x > y = case evalOp (x `Op.compare` y) :: GeneralDecimal of - Num { coefficient = 0 } -> False - Num { sign = Pos } -> True - _ -> False + x > y = case compareDecimal x y of + Right GT -> True + _ -> False - x >= y = case evalOp (x `Op.compare` y) :: GeneralDecimal of - Num { sign = Pos } -> 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) -{- $doctest-Ord -prop> x > y ==> max x y == x && max y x == (x :: BasicDecimal) -prop> x < y ==> min x y == x && min y x == (x :: BasicDecimal) - -prop> max x y == x ==> x >= y -prop> max x y == y ==> y >= x -prop> min x y == x ==> x <= y -prop> min x y == y ==> y <= x --} - -- | Unlike the instances for 'Float' and 'Double', the lists returned by the -- 'enumFromTo' and 'enumFromThenTo' methods in this instance terminate with -- the last element strictly less than (greater than in the case of a negative @@ -204,17 +213,6 @@ => Decimal p r -> Decimal p r -> [Decimal p r] enumFromWith x i = x : enumFromWith (x + i) i -{- $doctest-Enum ->>> [0, 0.1 .. 2] :: [BasicDecimal] -[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0] - ->>> [2, 1.9 .. 0] :: [BasicDecimal] -[2,1.9,1.8,1.7,1.6,1.5,1.4,1.3,1.2,1.1,1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1,0.0] - ->>> [1.7 .. 5.7] :: [BasicDecimal] -[1.7,2.7,3.7,4.7,5.7] --} - instance (Precision p, Rounding r) => Num (Decimal p r) where x + y = evalOp (x `Op.add` y) x - y = evalOp (x `Op.subtract` y) @@ -224,10 +222,11 @@ abs = evalOp . Op.abs signum n = case n of - Num { coefficient = 0 } -> zero - Num { sign = s } -> one { sign = s } - Inf { sign = s } -> one { sign = s } - _ -> n + 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 @@ -235,15 +234,6 @@ , exponent = 0 } -{- $doctest-Num -prop> x + x == x * (2 :: GeneralDecimal) -prop> isFinite x ==> x - x == (0 :: GeneralDecimal) -prop> isFinite x ==> x + negate x == (0 :: GeneralDecimal) -prop> abs x >= (0 :: GeneralDecimal) - -prop> abs x * signum x == (x :: GeneralDecimal) --} - 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) @@ -258,16 +248,6 @@ d = fromInteger (denominator r) :: GeneralDecimal in evalOp (n `Op.divide` d) -{- $doctest-Fractional -prop> (4.14 :: Decimal P2 RoundHalfUp) == 4.1 -prop> (4.15 :: Decimal P2 RoundHalfUp) == 4.2 -prop> (4.15 :: Decimal P2 RoundHalfDown) == 4.1 -prop> (4.15 :: Decimal P2 RoundHalfEven) == 4.2 -prop> (4.25 :: Decimal P2 RoundHalfEven) == 4.2 -prop> (4.35 :: Decimal P2 RoundHalfEven) == 4.4 -prop> (4.45 :: Decimal P2 RoundHalfEven) == 4.4 --} - instance (FinitePrecision p, Rounding r) => RealFrac (Decimal p r) where properFraction x@Num { sign = s, coefficient = c, exponent = e } | e < 0 = (n, f) @@ -277,95 +257,254 @@ (q, r) = c `quotRem` (10^(-e)) properFraction nan = (0, nan) -{- $doctest-RealFrac -prop> let (n,f) = properFraction (x :: BasicDecimal) in x == fromIntegral n + f -prop> let (n,f) = properFraction (x :: BasicDecimal) in (x < 0 && n <= 0) || (x >= 0 && n >= 0) -prop> let (n,f) = properFraction (x :: BasicDecimal) in (x < 0 && f <= 0) || (x >= 0 && f >= 0) -prop> let (n,f) = properFraction (x :: BasicDecimal) in isFinite f ==> abs f < 1 --} +-- | Compute a generalized continued fraction to maximum precision. A hint is +-- used to indicate the minimum number of terms that should be generated +-- before (expensively) examining the results for convergence. +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 + -- | Compute an infinite series to maximum precision. -infiniteSeries :: (FinitePrecision p, Rounding r) - => [Decimal p r] -> Decimal p r -infiniteSeries = series zero +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 + x + where n' = n `op` x series n [] = n --- | Compute the arcsine of the argument to maximum precision using series --- expansion. -arcsine :: (FinitePrecision p, Rounding r) => Decimal p r -> Decimal p r -arcsine x = infiniteSeries (x : series 1 2 x 3) - where series n d x i = - let x' = x * x2 - in (n * x') / (d * i) : series (n * i) (d * (i + one)) x' (i + two) +-- | Compute the inverse tangent of the argument to maximum precision using +-- series expansion. +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' + +-- | Compute the inverse sine of the argument to maximum precision using +-- series expansion. +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' + +-- | Compute the inverse tangent of the (decimal) argument to maximum +-- precision. +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 + +-- | Compute the inverse tangent of the (rational) argument to maximum +-- precision using a generalized continued fraction. +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 -- estimated minimum # terms + x2 = x * x + ty = 2 * y --- | Compute π to maximum precision using the arcsine series expansion. -seriesPi :: FinitePrecision p => Decimal p RoundHalfEven -seriesPi = 6 * arcsine oneHalf + -- [ (nx * nx, (n * 2 + 1) * y) | n <- [1..], let nx = n * x ] + partials n a b = + let a' = a + n * x2 + b' = b + ty + in (a', b') : partials (n + 2) a' b' +-- | Compute the inverse sine of the argument to maximum precision. +arcsin :: FinitePrecision p => Decimal p r -> ExtendedDecimal p +arcsin z = let z' = castUp z + in castDown' $ two * arctan (z' / (one + sqrt (one - z' * z'))) + +-- | Compute the inverse cosine of the argument to maximum precision. +arccos :: FinitePrecision p => Decimal p r -> ExtendedDecimal p +arccos = castDown' . (halfPi -) . arcsin . castUp + +-- | Compute π to maximum precision using the inverse sine series expansion. +seriesPi :: FinitePrecision p => ExtendedDecimal p +seriesPi = castDown' $ 6 * arcsin oneHalf + +-- | Compute π to maximum precision using the best-known Machin-like formula. +machinPi :: FinitePrecision p => ExtendedDecimal p +machinPi = castDown' $ 16 * arctan' (1 % 5) - 4 * arctan' (1 % 239) + +-- | Compute π to maximum precision using a generalized continued fraction +-- that converges linearly, adding at least three decimal digits of precision +-- per four terms. +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 + +-- | Precomputed π to a precision of 50 digits +fastPi :: FinitePrecision p => ExtendedDecimal p +fastPi = 3.1415926535897932384626433832795028841971693993751 + +-- | Compute π/2 to maximum precision. +halfPi :: FinitePrecision p => ExtendedDecimal p +halfPi = castDown' $ pi * oneHalf + +-- | Compute π/4 to maximum precision. +quarterPi :: FinitePrecision p => ExtendedDecimal p +quarterPi = castDown' $ pi * oneQuarter + +-- | Compute (cos 𝛽, sin 𝛽) to maximum precision using Volder's algorithm +-- (CORDIC). +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 = lim {n→∞} K(n) + -- K(n) = prod {i=0..n-1} 1 / sqrt (1 + 2^(-2 * i)) + 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) + +-- | Cast a number to a number with two additional digits of precision and +-- rounding half even. +castUp :: Precision p => Decimal p r -> ExtendedDecimal (PPlus1 (PPlus1 p)) +castUp = coerce + -- | Cast a number with two additional digits of precision down to a number --- with the desired precision. -castDown :: (Precision p, Rounding r) - => Decimal (PPlus1 (PPlus1 p)) a -> Decimal p r -castDown = cast +-- with the desired precision, rounding half even. +castDown' :: Precision p + => ExtendedDecimal (PPlus1 (PPlus1 p)) -> ExtendedDecimal p +castDown' = cast -notyet :: String -> a -notyet = error . (++ ": not yet implemented") +-- | Cast a number with two additional digits of precision down to a number +-- with the desired precision, rounding half even, but returning a number type +-- with arbitrary rounding. +castDown :: (Precision p, Rounding r) + => ExtendedDecimal (PPlus1 (PPlus1 p)) -> Decimal p r +castDown = castRounding . castDown' --- | The trigonometric and hyperbolic 'Floating' methods (other than the --- precision-dependent constant 'pi') are not yet implemented. +-- | The constant 'pi' is precision-dependent. instance (FinitePrecision p, Rounding r) => Floating (Decimal p r) where - pi = castDown seriesPi + 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 10 x = castRounding $ evalOp (Op.log10 x) - logBase _ 1 = zero - logBase b x = evalOp (join $ Op.divide <$> Op.ln x <*> Op.ln b) + 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 = notyet "sin" - cos = notyet "cos" - - asin = notyet "asin" - acos = notyet "acos" - atan = notyet "atan" - - sinh = notyet "sinh" - cosh = notyet "cosh" - - asinh = notyet "asinh" - acosh = notyet "acosh" - atanh = notyet "atanh" + sin = castDown . snd . cordic . castUp + cos = castDown . fst . cordic . castUp + tan = castDown . uncurry (flip (/)) . cordic . castUp -{- $doctest-Floating -prop> realToFrac (pi :: ExtendedDecimal P16) == (pi :: Double) + asin = castRounding . arcsin + acos = castRounding . arccos + atan = castRounding . arctan -prop> y >= 0 ==> (x :: BasicDecimal) ** fromInteger y == x ^ y + -- sinh x = let ex = exp x in (ex^2 - 1) / (2 * ex) + 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 = let ex = exp x in (ex^2 + 1) / (2 * ex) + 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 = let e2x = exp (2 * x) in (e2x - 1) / (e2x + 1) + tanh x = castDown . evalOp' $ + two `Op.multiply` x >>= Op.exp >>= \e2x -> + e2x `Op.subtract` one >>= \e2xm1 -> e2x `Op.add` one >>= (e2xm1 `Op.divide`) -prop> isFinite x && x >= 0 ==> coefficient (sqrt (x * x) - (x :: ExtendedDecimal P16)) <= 1 --} + -- asinh x = log (x + sqrt (x^2 + 1)) + asinh x = castDown . evalOp' $ x `Op.multiply` x >>= + (`Op.add` one) >>= Op.squareRoot >>= (x `Op.add`) >>= Op.ln + -- acosh x = log (x + sqrt (x^2 - 1)) + acosh x = castDown . evalOp' $ x `Op.multiply` x >>= + (`Op.subtract` one) >>= Op.squareRoot >>= (x `Op.add`) >>= Op.ln + -- atanh x = log ((1 + x) / (1 - x)) / 2 + 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 _ = (minBound, maxBound) -- ? + 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 } -> (m, n) - where m = signFunc s (fromIntegral c) - n = fromIntegral e - Inf { sign = s } -> (special s 0, maxBound ) - QNaN { sign = s, payload = p } -> (special s p, minBound ) - SNaN { sign = s, payload = p } -> (special s p, minBound + 1) + 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 @@ -378,16 +517,16 @@ } special | n == maxBound = Inf { sign = signMatch m } - | n == minBound = QNaN { sign = signMatch m, payload = p } - | otherwise = SNaN { sign = signMatch m, payload = p } + | 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 - QNaN{} -> True - SNaN{} -> True - _ -> False + NaN{} -> True + _ -> False isInfinite x = case x of Inf{} -> True @@ -401,17 +540,6 @@ isIEEE _ = True -{- $doctest-RealFloat -prop> uncurry encodeFloat (decodeFloat x) == (x :: BasicDecimal) -prop> isFinite x ==> significand x * fromInteger (floatRadix x) ^^ Prelude.exponent x == (x :: BasicDecimal) -prop> Prelude.exponent (0 :: BasicDecimal) == 0 -prop> isFinite x && x /= 0 ==> Prelude.exponent (x :: BasicDecimal) == snd (decodeFloat x) + floatDigits x - -prop> isNegativeZero (read "-0" :: BasicDecimal) == True -prop> isNegativeZero (read "+0" :: BasicDecimal) == False -prop> x /= 0 ==> isNegativeZero (x :: BasicDecimal) == False --} - -- | The 'Bits' instance makes use of the logical operations from the -- /General Decimal Arithmetic Specification/ using a /digit-wise/ -- representation of bits where the /sign/ is non-negative, the /exponent/ is @@ -454,6 +582,13 @@ 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 + -- | A 'Decimal' representing the value zero zero :: Decimal p r zero = Num { sign = Pos @@ -465,6 +600,10 @@ oneHalf :: Decimal p r oneHalf = zero { coefficient = 5, exponent = -1 } +-- | A 'Decimal' representing the value ¼ +oneQuarter :: Decimal p r +oneQuarter = zero { coefficient = 25, exponent = -2 } + -- | A 'Decimal' representing the value one one :: Decimal p r one = zero { coefficient = 1 } @@ -473,6 +612,10 @@ two :: Decimal p r two = zero { coefficient = 2 } +-- | A 'Decimal' representing the value three +three :: Decimal p r +three = zero { coefficient = 3 } + -- | A 'Decimal' representing the value ten ten :: Decimal p r ten = zero { coefficient = 10 } @@ -487,11 +630,11 @@ -- | A 'Decimal' representing undefined results qNaN :: Decimal p r -qNaN = QNaN { sign = Pos, payload = 0 } +qNaN = NaN { sign = Pos, signaling = False, payload = 0 } -- | A signaling 'Decimal' representing undefined results sNaN :: Decimal p r -sNaN = SNaN { sign = Pos, payload = 0 } +sNaN = qNaN { signaling = True } -- | Negate the given 'Decimal' by directly flipping its sign. flipSign :: Decimal p r -> Decimal p r @@ -559,35 +702,31 @@ | isFinite n && not (isZero n) = maybe False (adjustedExponent n <) (eMin n) | otherwise = False --- | If the argument is 'False', return a 'Decimal' value zero; if 'True', --- return the value one. This is basically an optimized @toEnum . fromEnum@ to --- support an all-decimal usage of the operations from --- "Numeric.Decimal.Operation" that return a 'Bool'. +-- | Return @0@ or @1@ if the argument is 'False' or 'True', respectively. +-- This is basically an optimized @'toEnum' . 'fromEnum'@ and allows an +-- all-decimal usage of the operations from "Numeric.Decimal.Operation" that +-- return a 'Bool'. fromBool :: Bool -> Decimal p r fromBool False = zero fromBool True = one -- | Return 'False' if the argument is zero or NaN, and 'True' otherwise. toBool :: Decimal p r -> Bool -toBool Num { coefficient = c } - | c == 0 = False - | otherwise = True -toBool Inf{} = True -toBool _ = False +toBool Num { coefficient = c } = c /= 0 +toBool NaN{} = False +toBool _ = True +-- | Return @-1@, @0@, or @1@ if the argument is 'LT', 'EQ', or 'GT', +-- respectively. This allows an all-decimal usage of the operations from +-- "Numeric.Decimal.Operation" that return an 'Ordering'. +fromOrdering :: Ordering -> Decimal p r +fromOrdering LT = negativeOne +fromOrdering EQ = zero +fromOrdering GT = one + -- | Upper limit on the absolute value of the exponent eLimit :: Precision p => p -> Maybe Exponent eLimit = eMax -- ? - --- | Minimum value of the adjusted exponent -eMin :: Precision p => p -> Maybe Exponent -eMin n = (1 -) <$> eMax n - --- | Maximum value of the adjusted exponent -eMax :: Precision p => p -> Maybe Exponent -eMax n = subtract 1 . (10 ^) . numDigits <$> base - where mlength = precision n :: Maybe Int - base = (10 *) . fromIntegral <$> mlength :: Maybe Coefficient -- | Minimum value of the exponent for subnormal results eTiny :: Precision p => p -> Maybe Exponent diff --git a/src/Numeric/Decimal/Number.hs-boot b/src/Numeric/Decimal/Number.hs-boot --- a/src/Numeric/Decimal/Number.hs-boot +++ b/src/Numeric/Decimal/Number.hs-boot @@ -4,37 +4,31 @@ module Numeric.Decimal.Number ( Sign(..) - , Decimal(..) , Coefficient + , Exponent + , Decimal(..) , numDigits ) where import Numeric.Natural (Natural) -import Numeric.Decimal.Precision - data Sign = Pos | Neg -instance Eq Sign type Coefficient = Natural -type Exponent = Int +type Exponent = Integer type Payload = Coefficient type role Decimal phantom phantom data Decimal p r - = Num { sign :: Sign - , coefficient :: Coefficient - , exponent :: Exponent - } - | Inf { sign :: Sign - } - | QNaN { sign :: Sign - , payload :: Payload - } - | SNaN { sign :: Sign - , payload :: Payload - } - -instance Precision p => Precision (Decimal p r) + = Num { sign :: Sign + , coefficient :: Coefficient + , exponent :: Exponent + } + | Inf { sign :: Sign + } + | NaN { sign :: Sign + , signaling :: Bool + , payload :: Payload + } numDigits :: Coefficient -> Int diff --git a/src/Numeric/Decimal/Operation.hs b/src/Numeric/Decimal/Operation.hs --- a/src/Numeric/Decimal/Operation.hs +++ b/src/Numeric/Decimal/Operation.hs @@ -16,1950 +16,1383 @@ ( -- * Arithmetic operations -- $arithmetic-operations - abs - , add - , subtract - , compare - , compareSignal - , divide - -- divideInteger - , exp - , fusedMultiplyAdd - , ln - , log10 - , max - , maxMagnitude - , min - , minMagnitude - , minus - , plus - , multiply - -- nextMinus - -- nextPlus - -- nextToward - , power - , quantize - , reduce - -- remainder - -- remainderNear - -- roundToIntegralExact - -- roundToIntegralValue - , squareRoot - - -- * Miscellaneous operations - -- $miscellaneous-operations - - , and - , canonical - , class_, Class(..), Sign(..), NumberClass(..), NaNClass(..) - -- compareTotal - -- compareTotalMagnitude - , copy - , copyAbs - , copyNegate - , copySign - , invert - , isCanonical - , isFinite - , isInfinite - , isNaN - , isNormal - , isQNaN - , isSigned - , isSNaN - , isSubnormal - , isZero - , logb - , or - , radix - , rotate - , sameQuantum - -- scaleb - , shift - , xor - ) where - -import Prelude hiding (abs, and, compare, exp, exponent, isInfinite, isNaN, - max, min, or, subtract) -import qualified Prelude - -import Control.Monad (join) -import Data.Bits (complement, setBit, testBit, zeroBits, (.&.), (.|.)) -import Data.Coerce (coerce) -import Data.List (find) -import Data.Maybe (fromMaybe) - -import qualified Data.Bits as Bits - -import Numeric.Decimal.Arithmetic -import Numeric.Decimal.Number hiding (isFinite, isNormal, isSubnormal, isZero) -import Numeric.Decimal.Precision -import Numeric.Decimal.Rounding - -import qualified Numeric.Decimal.Number as Number - -{- $setup ->>> :load Harness --} - -finitePrecision :: FinitePrecision p => Decimal p r -> Int -finitePrecision n = let Just p = precision n in p - -roundingAlg :: Rounding r => Arith p r a -> RoundingAlgorithm -roundingAlg = rounding . arithRounding - where arithRounding :: Arith p r a -> r - arithRounding = undefined - -result :: (Precision p, Rounding r) => Decimal p r -> Arith p r (Decimal p r) -result = roundDecimal -- ... --- | maybe False (numDigits c >) (precision r) = undefined - -invalidOperation :: Decimal a b -> Arith p r (Decimal p r) -invalidOperation n = raiseSignal InvalidOperation qNaN - -toQNaN :: Decimal a b -> Decimal p r -toQNaN SNaN { sign = s, payload = p } = QNaN { sign = s, payload = p } -toQNaN n@QNaN{} = coerce n -toQNaN n = qNaN { sign = sign n } - -toQNaN2 :: Decimal a b -> Decimal c d -> Decimal p r -toQNaN2 nan@SNaN{} _ = toQNaN nan -toQNaN2 _ nan@SNaN{} = toQNaN nan -toQNaN2 nan@QNaN{} _ = coerce nan -toQNaN2 _ nan@QNaN{} = coerce nan -toQNaN2 n _ = toQNaN n - -quietToSignal :: Decimal p r -> Decimal p r -quietToSignal QNaN { sign = s, payload = p } = SNaN { sign = s, payload = p } -quietToSignal x = x - --- $arithmetic-operations --- --- This section describes the arithmetic operations on, and some other --- functions of, numbers, including subnormal numbers, negative zeros, and --- special values (see also IEEE 754 §5 and §6). - -{- $doctest-special-values ->>> op2 Op.add "Infinity" "1" -Infinity - ->>> op2 Op.add "NaN" "1" -NaN - ->>> op2 Op.add "NaN" "Infinity" -NaN - ->>> op2 Op.subtract "1" "Infinity" --Infinity - ->>> op2 Op.multiply "-1" "Infinity" --Infinity - ->>> op2 Op.subtract "-0" "0" --0 - ->>> op2 Op.multiply "-1" "0" --0 - ->>> op2 Op.divide "1" "0" -Infinity - ->>> op2 Op.divide "1" "-0" --Infinity - ->>> op2 Op.divide "-1" "0" --Infinity --} - --- | 'add' takes two operands. If either operand is a /special value/ then the --- general rules apply. --- --- Otherwise, the operands are added. --- --- The result is then rounded to /precision/ digits if necessary, counting --- from the most significant digit of the result. -add :: (Precision p, Rounding r) - => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) -add Num { sign = xs, coefficient = xc, exponent = xe } - Num { sign = ys, coefficient = yc, exponent = ye } = sum - - where sum = result Num { sign = rs, coefficient = rc, exponent = re } - rs | rc /= 0 = if xac > yac then xs else ys - | xs == Neg && ys == Neg = Neg - | xs /= ys && - roundingAlg sum == RoundFloor = Neg - | otherwise = Pos - rc | xs == ys = xac + yac - | xac > yac = xac - yac - | otherwise = yac - xac - re = Prelude.min xe ye - (xac, yac) | xe == ye = (xc, yc) - | xe > ye = (xc * 10^n, yc) - | otherwise = (xc, yc * 10^n) - where n = Prelude.abs (xe - ye) - -add inf@Inf { sign = xs } Inf { sign = ys } - | xs == ys = return (coerce inf) - | otherwise = invalidOperation inf -add inf@Inf{} Num{} = return (coerce inf) -add Num{} inf@Inf{} = return (coerce inf) -add x y = return (toQNaN2 x y) - -{- $doctest-add ->>> op2 Op.add "12" "7.00" -19.00 - ->>> op2 Op.add "1E+2" "1E+4" -1.01E+4 --} - --- | 'subtract' takes two operands. If either operand is a /special value/ --- then the general rules apply. --- --- Otherwise, the operands are added after inverting the /sign/ used for the --- second operand. --- --- The result is then rounded to /precision/ digits if necessary, counting --- from the most significant digit of the result. -subtract :: (Precision p, Rounding r) - => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) -subtract x = add x . flipSign - -{- $doctest-subtract ->>> op2 Op.subtract "1.3" "1.07" -0.23 - ->>> op2 Op.subtract "1.3" "1.30" -0.00 - ->>> op2 Op.subtract "1.3" "2.07" --0.77 --} - --- | 'minus' takes one operand, and corresponds to the prefix minus operator --- in programming languages. --- --- Note that the result of this operation is affected by context and may set --- /flags/. The 'copyNegate' operation may be used instead of 'minus' if this --- is not desired. -minus :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r) -minus x = zero { exponent = exponent x } `subtract` x - -{- $doctest-minus ->>> op1 Op.minus "1.3" --1.3 - ->>> op1 Op.minus "-1.3" -1.3 --} - --- | 'plus' takes one operand, and corresponds to the prefix plus operator in --- programming languages. --- --- Note that the result of this operation is affected by context and may set --- /flags/. -plus :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r) -plus x = zero { exponent = exponent x } `add` x - -{- $doctest-plus ->>> op1 Op.plus "1.3" -1.3 - ->>> op1 Op.plus "-1.3" --1.3 --} - --- | 'multiply' takes two operands. If either operand is a /special value/ --- then the general rules apply. Otherwise, the operands are multiplied --- together (“long multiplication”), resulting in a number which may be as --- long as the sum of the lengths of the two operands. --- --- The result is then rounded to /precision/ digits if necessary, counting --- from the most significant digit of the result. -multiply :: (Precision p, Rounding r) - => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) -multiply Num { sign = xs, coefficient = xc, exponent = xe } - Num { sign = ys, coefficient = yc, exponent = ye } = result rn - - where rn = Num { sign = rs, coefficient = rc, exponent = re } - rs = xorSigns xs ys - rc = xc * yc - re = xe + ye - -multiply Inf { sign = xs } Inf { sign = ys } = - return Inf { sign = xorSigns xs ys } -multiply Inf { sign = xs } Num { sign = ys, coefficient = yc } - | yc == 0 = invalidOperation qNaN - | otherwise = return Inf { sign = xorSigns xs ys } -multiply Num { sign = xs, coefficient = xc } Inf { sign = ys } - | xc == 0 = invalidOperation qNaN - | otherwise = return Inf { sign = xorSigns xs ys } -multiply nan@SNaN{} _ = invalidOperation nan -multiply _ nan@SNaN{} = invalidOperation nan -multiply x y = return (toQNaN2 x y) - -{- $doctest-multiply ->>> op2 Op.multiply "1.20" "3" -3.60 - ->>> op2 Op.multiply "7" "3" -21 - ->>> op2 Op.multiply "0.9" "0.8" -0.72 - ->>> op2 Op.multiply "0.9" "-0" --0.0 - ->>> op2 Op.multiply "654321" "654321" -4.28135971E+11 --} - --- | 'exp' takes one operand. If the operand is a NaN then the general rules --- for special values apply. --- --- Otherwise, the result is /e/ raised to the power of the operand, with the --- following cases: --- --- * If the operand is −Infinity, the result is 0 and exact. --- --- * If the operand is a zero, the result is 1 and exact. --- --- * If the operand is +Infinity, the result is +Infinity and exact. --- --- * Otherwise the result is inexact and will be rounded using the --- /round-half-even/ algorithm. The coefficient will have exactly /precision/ --- digits (unless the result is subnormal). These inexact results should be --- correctly rounded, but may be up to 1 ulp (unit in last place) in error. -exp :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p RoundHalfEven) -exp x@Num { sign = s, coefficient = c } - | c == 0 = return one - | s == Neg = subArith (maclaurin x { sign = Pos } >>= reciprocal) >>= - subRounded >>= result - | otherwise = subArith (maclaurin x) >>= subRounded >>= result - - where multiplyExact :: Decimal a b -> Decimal c d - -> Arith PInfinite RoundHalfEven - (Decimal PInfinite RoundHalfEven) - multiplyExact = multiply - - maclaurin :: FinitePrecision p => Decimal a b - -> Arith p RoundHalfEven (Decimal p RoundHalfEven) - maclaurin x - | adjustedExponent x >= 0 = subArith (subMaclaurin x) >>= subRounded - | otherwise = sum one one one one - where sum :: FinitePrecision p - => Decimal p RoundHalfEven - -> Decimal PInfinite RoundHalfEven - -> Decimal PInfinite RoundHalfEven - -> Decimal PInfinite RoundHalfEven - -> Arith p RoundHalfEven (Decimal p RoundHalfEven) - sum s num den n = do - num' <- subArith (multiplyExact num x) - den' <- subArith (multiplyExact den n) - s' <- add s =<< divide num' den' - if s' == s then return s' - else sum s' num' den' =<< subArith (add n one) - - subMaclaurin :: FinitePrecision p => Decimal a b - -> Arith p RoundHalfEven (Decimal p RoundHalfEven) - subMaclaurin x = subArith (multiplyExact x oneHalf) >>= maclaurin >>= - \r -> multiply r r - - subRounded :: Precision p - => Decimal (PPlus1 (PPlus1 p)) a - -> Arith p r (Decimal p RoundHalfEven) - subRounded = subArith . roundDecimal - - result :: Decimal p a -> Arith p r (Decimal p a) - result r = coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r') - where r' = coerce r - -exp n@Inf { sign = s } - | s == Pos = return (coerce n) - | otherwise = return zero -exp n@QNaN{} = return (coerce n) -exp n@SNaN{} = coerce <$> invalidOperation n - -{- $doctest-exp ->>> op1 Op.exp "-Infinity" -0 - ->>> op1 Op.exp "-1" -0.367879441 - ->>> op1 Op.exp "0" -1 - ->>> op1 Op.exp "1" -2.71828183 - ->>> op1 Op.exp "0.693147181" -2.00000000 - ->>> op1 Op.exp "+Infinity" -Infinity --} - --- | 'fusedMultiplyAdd' takes three operands; the first two are multiplied --- together, using 'multiply', with sufficient precision and exponent range --- that the result is exact and unrounded. No /flags/ are set by the --- multiplication unless one of the first two operands is a signaling NaN or --- one is a zero and the other is an infinity. --- --- Unless the multiplication failed, the third operand is then added to the --- result of that multiplication, using 'add', under the current context. --- --- In other words, @fusedMultiplyAdd x y z@ delivers a result which is @(x × --- y) + z@ with only the one, final, rounding. -fusedMultiplyAdd :: (Precision p, Rounding r) - => Decimal a b -> Decimal c d -> Decimal e f - -> Arith p r (Decimal p r) -fusedMultiplyAdd x y z = - either raise (return . coerce) (exactMult x y) >>= add z - - where exactMult :: Rounding r => Decimal a b -> Decimal c d - -> Either (Exception PInfinite r) (Decimal PInfinite r) - exactMult x y = evalArith (multiply x y) newContext - - raise :: Exception a r -> Arith p r (Decimal p r) - raise e = raiseSignal (exceptionSignal e) (coerce $ exceptionResult e) - -{- $doctest-fusedMultiplyAdd ->>> op3 Op.fusedMultiplyAdd "3" "5" "7" -22 - ->>> op3 Op.fusedMultiplyAdd "3" "-5" "7" --8 - ->>> op3 Op.fusedMultiplyAdd "888565290" "1557.96930" "-86087.7578" -1.38435736E+12 --} - --- | 'ln' takes one operand. If the operand is a NaN then the general rules --- for special values apply. --- --- Otherwise, the operand must be a zero or positive, and the result is the --- natural (base /e/) logarithm of the operand, with the following cases: --- --- * If the operand is a zero, the result is −Infinity and exact. --- --- * If the operand is +Infinity, the result is +Infinity and exact. --- --- * If the operand equals one, the result is 0 and exact. --- --- * Otherwise the result is inexact and will be rounded using the --- /round-half-even/ algorithm. The coefficient will have exactly /precision/ --- digits (unless the result is subnormal). These inexact results should be --- correctly rounded, but may be up to 1 ulp (unit in last place) in error. -ln :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p RoundHalfEven) -ln x@Num { sign = s, coefficient = c, exponent = e } - | c == 0 = return infinity { sign = Neg } - | s == Pos = if e <= 0 && c == 10^(-e) then return zero - else subArith (subLn x) >>= subRounded >>= result - - where subLn :: FinitePrecision p => Decimal a b - -> Arith p RoundHalfEven (Decimal p RoundHalfEven) - subLn x = do - let fe = fromIntegral (-(numDigits c - 1)) :: Exponent - r = fromIntegral (e - fe) :: Decimal PInfinite RoundHalfEven - lnf <- taylorLn x { exponent = fe } - add lnf =<< multiply r =<< ln10 - - subRounded :: Precision p => Decimal (PPlus1 (PPlus1 p)) a - -> Arith p r (Decimal p RoundHalfEven) - subRounded = subArith . roundDecimal - - result :: Decimal p a -> Arith p r (Decimal p a) - result r = coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r') - where r' = coerce r - -ln n@Inf { sign = Pos } = return (coerce n) -ln n@QNaN{} = return (coerce n) -ln n = coerce <$> invalidOperation n - -taylorLn :: FinitePrecision p => Decimal a b - -> Arith p RoundHalfEven (Decimal p RoundHalfEven) -taylorLn x = do - num <- x `subtract` one - den <- x `add` one - multiply two =<< sum =<< num `divide` den - - where sum :: FinitePrecision p => Decimal p RoundHalfEven - -> Arith p RoundHalfEven (Decimal p RoundHalfEven) - sum b = multiply b b >>= \b2 -> sum' b b b2 one - - where sum' :: FinitePrecision p - => Decimal p RoundHalfEven - -> Decimal p RoundHalfEven - -> Decimal p RoundHalfEven - -> Decimal PInfinite RoundHalfEven - -> Arith p RoundHalfEven (Decimal p RoundHalfEven) - sum' s m b n = do - m' <- multiply m b - n' <- subArith (add n two) - s' <- add s =<< divide m' n' - if s' == s then return s' else sum' s' m' b n' - -ln10 :: FinitePrecision p => Arith p r (Decimal p RoundHalfEven) -ln10 = getPrecision >>= \(Just p) -> - if p <= 50 then return fastLn10 else slowLn10 - - where fastLn10 :: FinitePrecision p => Decimal p RoundHalfEven - fastLn10 = 2.3025850929940456840179914546843642076011014886288 - - slowLn10 :: FinitePrecision p => Arith p r (Decimal p RoundHalfEven) - slowLn10 = subArith (taylorLn ten) >>= subRound - - where subRound :: Precision p => Decimal (PPlus1 (PPlus1 p)) a - -> Arith p r (Decimal p RoundHalfEven) - subRound = subArith . roundDecimal - -{- $doctest-ln ->>> op1 Op.ln "0" --Infinity - ->>> op1 Op.ln "1.000" -0 - ->>> op1 Op.ln "2.71828183" -1.00000000 - ->>> op1 Op.ln "10" -2.30258509 - ->>> op1 Op.ln "+Infinity" -Infinity --} - --- | 'log10' takes one operand. If the operand is a NaN then the general rules --- for special values apply. --- --- Otherwise, the operand must be a zero or positive, and the result is the --- base 10 logarithm of the operand, with the following cases: --- --- * If the operand is a zero, the result is −Infinity and exact. --- --- * If the operand is +Infinity, the result is +Infinity and exact. --- --- * If the operand equals an integral power of ten (including 10^0 and --- negative powers) and there is sufficient /precision/ to hold the integral --- part of the result, the result is an integer (with an exponent of 0) and --- exact. --- --- * Otherwise the result is inexact and will be rounded using the --- /round-half-even/ algorithm. The coefficient will have exactly /precision/ --- digits (unless the result is subnormal). These inexact results should be --- correctly rounded, but may be up to 1 ulp (unit in last place) in error. -log10 :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p RoundHalfEven) -log10 x@Num { sign = s, coefficient = c, exponent = e } - | c == 0 = return infinity { sign = Neg } - | s == Pos = getPrecision >>= \prec -> case powerOfTen c e of - Just p | maybe True (numDigits pc <=) prec -> return (fromInteger p) - where pc = fromInteger (Prelude.abs p) :: Coefficient - _ -> subArith (join $ divide <$> ln x <*> ln10) >>= result - - where powerOfTen :: Coefficient -> Exponent -> Maybe Integer - powerOfTen c e - | c == 10^d = Just (fromIntegral e + fromIntegral d) - | otherwise = Nothing - where d = numDigits c - 1 :: Int - - result :: Decimal p a -> Arith p r (Decimal p a) - result r = coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r') - where r' = coerce r - -log10 n@Inf { sign = Pos } = return (coerce n) -log10 n@QNaN{} = return (coerce n) -log10 n = coerce <$> invalidOperation n - -{- $doctest-log10 ->>> op1 Op.log10 "0" --Infinity - ->>> op1 Op.log10 "0.001" --3 - ->>> op1 Op.log10 "1.000" -0 - ->>> op1 Op.log10 "2" -0.301029996 - ->>> op1 Op.log10 "10" -1 - ->>> op1 Op.log10 "70" -1.84509804 - ->>> op1 Op.log10 "+Infinity" -Infinity --} - --- | 'divide' takes two operands. If either operand is a /special value/ then --- the general rules apply. --- --- Otherwise, if the divisor is zero then either the Division undefined --- condition is raised (if the dividend is zero) and the result is NaN, or the --- Division by zero condition is raised and the result is an Infinity with a --- sign which is the exclusive or of the signs of the operands. --- --- Otherwise, a “long division” is effected. --- --- The result is then rounded to /precision/ digits, if necessary, according --- to the /rounding/ algorithm and taking into account the remainder from the --- division. -divide :: (FinitePrecision p, Rounding r) - => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) -divide dividend@Num{ sign = xs } Num { coefficient = 0, sign = ys } - | Number.isZero dividend = invalidOperation qNaN - | otherwise = raiseSignal DivisionByZero - infinity { sign = xorSigns xs ys } -divide Num { sign = xs, coefficient = xc, exponent = xe } - Num { sign = ys, coefficient = yc, exponent = ye } = quotient - - where quotient = result =<< answer - rn = Num { sign = rs, coefficient = rc, exponent = re } - rs = xorSigns xs ys - (rc, rem, dv, adjust) = longDivision xc yc (finitePrecision rn) - re = xe - (ye + adjust) - answer - | rem == 0 = return rn - | otherwise = roundDecimal $ case (rem * 2) `Prelude.compare` dv of - LT -> rn { coefficient = rc * 10 + 1, exponent = re - 1 } - EQ -> rn { coefficient = rc * 10 + 5, exponent = re - 1 } - GT -> rn { coefficient = rc * 10 + 9, exponent = re - 1 } - -divide Inf{} Inf{} = invalidOperation qNaN -divide Inf { sign = xs } Num { sign = ys } = - return Inf { sign = xorSigns xs ys } -divide Num { sign = xs } Inf { sign = ys } = - return zero { sign = xorSigns xs ys } -divide x y = return (toQNaN2 x y) - -{- $doctest-divide ->>> op2 Op.divide "1" "3" -0.333333333 - ->>> op2 Op.divide "2" "3" -0.666666667 - ->>> op2 Op.divide "5" "2" -2.5 - ->>> op2 Op.divide "1" "10" -0.1 - ->>> op2 Op.divide "12" "12" -1 - ->>> op2 Op.divide "8.00" "2" -4.00 - ->>> op2 Op.divide "2.400" "2.0" -1.20 - ->>> op2 Op.divide "1000" "100" -10 - ->>> op2 Op.divide "1000" "1" -1000 - ->>> op2 Op.divide "2.40E+6" "2" -1.20E+6 --} - -type Dividend = Coefficient -type Divisor = Coefficient -type Quotient = Coefficient -type Remainder = Dividend - -longDivision :: Dividend -> Divisor -> Int - -> (Quotient, Remainder, Divisor, Exponent) -longDivision 0 dv _ = (0, 0, dv, 0) -longDivision dd dv p = step1 dd dv 0 - - where step1 :: Dividend -> Divisor -> Exponent - -> (Quotient, Remainder, Divisor, Exponent) - step1 dd dv adjust - | dd < dv = step1 (dd * 10) dv (adjust + 1) - | dd >= 10 * dv = step1 dd (dv * 10) (adjust - 1) - | otherwise = step2 dd dv adjust - - step2 :: Dividend -> Divisor -> Exponent - -> (Quotient, Remainder, Divisor, Exponent) - step2 = step3 0 - - step3 :: Quotient -> Dividend -> Divisor -> Exponent - -> (Quotient, Remainder, Divisor, Exponent) - step3 r dd dv adjust - | dv <= dd = step3 (r + 1) (dd - dv) dv adjust - | (dd == 0 && adjust >= 0) || - numDigits r == p = step4 r dd dv adjust - | otherwise = step3 (r * 10) (dd * 10) dv (adjust + 1) - - step4 :: Quotient -> Remainder -> Divisor -> Exponent - -> (Quotient, Remainder, Divisor, Exponent) - step4 = (,,,) - -reciprocal :: (FinitePrecision p, Rounding r) - => Decimal a b -> Arith p r (Decimal p r) -reciprocal = divide one - --- | 'abs' takes one operand. If the operand is negative, the result is the --- same as using the 'minus' operation on the operand. Otherwise, the result --- is the same as using the 'plus' operation on the operand. --- --- Note that the result of this operation is affected by context and may set --- /flags/. The 'copyAbs' operation may be used if this is not desired. -abs :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r) -abs x - | isNegative x = minus x - | otherwise = plus x - -{- $doctest-abs ->>> op1 Op.abs "2.1" -2.1 - ->>> op1 Op.abs "-100" -100 - ->>> op1 Op.abs "101.5" -101.5 - ->>> op1 Op.abs "-101.5" -101.5 --} - --- | 'compare' takes two operands and compares their values numerically. If --- either operand is a /special value/ then the general rules apply. No flags --- are set unless an operand is a signaling NaN. --- --- Otherwise, the operands are compared, returning @−1@ if the first is less --- than the second, @0@ if they are equal, or @1@ if the first is greater than --- the second. -compare :: (Precision p, Rounding r) - => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) -compare x@Num{} y@Num{} = nzp <$> (xn `subtract` yn) - - where (xn, yn) | sign x /= sign y = (nzp x, nzp y) - | otherwise = (x, y) - - nzp :: Decimal p r -> Decimal p r - nzp Num { sign = s, coefficient = c } - | c == 0 = zero - | s == Pos = one - | otherwise = negativeOne - nzp Inf { sign = s } - | s == Pos = one - | otherwise = negativeOne - nzp n = toQNaN n - -compare Inf { sign = xs } Inf { sign = ys } - | xs == ys = return zero - | xs == Neg = return negativeOne - | otherwise = return one -compare Inf { sign = xs } Num { } - | xs == Neg = return negativeOne - | otherwise = return one -compare Num { } Inf { sign = ys } - | ys == Pos = return negativeOne - | otherwise = return one -compare nan@SNaN{} _ = invalidOperation nan -compare _ nan@SNaN{} = invalidOperation nan -compare x y = return (toQNaN2 x y) - -{- $doctest-compare ->>> op2 Op.compare "2.1" "3" --1 - ->>> op2 Op.compare "2.1" "2.1" -0 - ->>> op2 Op.compare "2.1" "2.10" -0 - ->>> op2 Op.compare "3" "2.1" -1 - ->>> op2 Op.compare "2.1" "-3" -1 - ->>> op2 Op.compare "-3" "2.1" --1 --} - --- | 'compareSignal' takes two operands and compares their values --- numerically. This operation is identical to 'compare', except that if --- neither operand is a signaling NaN then any quiet NaN operand is treated as --- though it were a signaling NaN. (That is, all NaNs signal, with signaling --- NaNs taking precedence over quiet NaNs.) -compareSignal :: (Precision p, Rounding r) - => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) -compareSignal x@SNaN{} y = x `compare` y -compareSignal x y@SNaN{} = x `compare` y -compareSignal x y = quietToSignal x `compare` quietToSignal y - --- | 'max' takes two operands, compares their values numerically, and returns --- the maximum. If either operand is a NaN then the general rules apply, --- unless one is a quiet NaN and the other is numeric, in which case the --- numeric operand is returned. -max :: (Precision p, Rounding r) - => Decimal a b -> Decimal a b -> Arith p r (Decimal a b) -max x y = snd <$> minMax id x y - -{- $doctest-max ->>> op2 Op.max "3" "2" -3 - ->>> op2 Op.max "-10" "3" -3 - ->>> op2 Op.max "1.0" "1" -1 - ->>> op2 Op.max "7" "NaN" -7 --} - --- | 'maxMagnitude' takes two operands and compares their values numerically --- with their /sign/ ignored and assumed to be 0. --- --- If, without signs, the first operand is the larger then the original first --- operand is returned (that is, with the original sign). If, without signs, --- the second operand is the larger then the original second operand is --- returned. Otherwise the result is the same as from the 'max' operation. -maxMagnitude :: (Precision p, Rounding r) - => Decimal a b -> Decimal a b -> Arith p r (Decimal a b) -maxMagnitude x y = snd <$> minMax withoutSign x y - --- | 'min' takes two operands, compares their values numerically, and returns --- the minimum. If either operand is a NaN then the general rules apply, --- unless one is a quiet NaN and the other is numeric, in which case the --- numeric operand is returned. -min :: (Precision p, Rounding r) - => Decimal a b -> Decimal a b -> Arith p r (Decimal a b) -min x y = fst <$> minMax id x y - -{- $doctest-min ->>> op2 Op.min "3" "2" -2 - ->>> op2 Op.min "-10" "3" --10 - ->>> op2 Op.min "1.0" "1" -1.0 - ->>> op2 Op.min "7" "NaN" -7 --} - --- | 'minMagnitude' takes two operands and compares their values numerically --- with their /sign/ ignored and assumed to be 0. --- --- If, without signs, the first operand is the smaller then the original first --- operand is returned (that is, with the original sign). If, without signs, --- the second operand is the smaller then the original second operand is --- returned. Otherwise the result is the same as from the 'min' operation. -minMagnitude :: (Precision p, Rounding r) - => Decimal a b -> Decimal a b -> Arith p r (Decimal a b) -minMagnitude x y = fst <$> minMax withoutSign x y - --- | Ordering function for 'min', 'minMagnitude', 'max', and 'maxMagnitude': --- returns the original arguments as (smaller, larger) when the given function --- is applied to them. -minMax :: (Precision p, Rounding r) - => (Decimal a b -> Decimal a b) -> Decimal a b -> Decimal a b - -> Arith p r (Decimal a b, Decimal a b) -minMax _ x@Num{} QNaN{} = return (x, x) -minMax _ x@Inf{} QNaN{} = return (x, x) -minMax _ QNaN{} y@Num{} = return (y, y) -minMax _ QNaN{} y@Inf{} = return (y, y) - -minMax f x y = do - c <- f x `compare` f y - return $ case c of - Num { coefficient = 0 } -> case (sign x, sign y) of - (Neg, Pos) -> (x, y) - (Pos, Neg) -> (y, x) - (Pos, Pos) -> case (x, y) of - (Num { exponent = xe }, Num { exponent = ye }) | xe > ye -> (y, x) - _ -> (x, y) - (Neg, Neg) -> case (x, y) of - (Num { exponent = xe }, Num { exponent = ye }) | xe < ye -> (y, x) - _ -> (x, y) - Num { sign = Pos } -> (y, x) - Num { sign = Neg } -> (x, y) - nan -> let nan' = coerce nan in (nan', nan') - -withoutSign :: Decimal p r -> Decimal p r -withoutSign n = n { sign = Pos } - --- | 'power' takes two operands, and raises a number (the left-hand operand) --- to a power (the right-hand operand). If either operand is a /special value/ --- then the general rules apply, except in certain cases. -power :: (FinitePrecision p, Rounding r) - => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) -power x@Num { coefficient = 0 } y@Num{} - | Number.isZero y = invalidOperation qNaN - | Number.isNegative y = return infinity { sign = powerSign x y } - | otherwise = return zero { sign = powerSign x y } -power x@Num{} y@Num{} = case integralValue y of - Just i | i < 0 -> reciprocal x >>= \rx -> integralPower rx (-i) - | otherwise -> integralPower x i - Nothing | Number.isPositive x -> ln x >>= multiply y >>= fmap coerce . exp - | otherwise -> invalidOperation qNaN -power x@Num{} y@Inf{} - | Number.isPositive x = return $ case sign y of - Pos -> infinity - Neg -> zero - | otherwise = invalidOperation qNaN -power x@Inf{} y@Num{} - | Number.isZero y = return one - | Number.isPositive y = return infinity { sign = powerSign x y } - | otherwise = return zero { sign = powerSign x y } -power Inf{} Inf { sign = s } - | s == Pos = return infinity - | otherwise = return zero -power x@SNaN{} _ = invalidOperation x -power _ y@SNaN{} = invalidOperation y -power x@QNaN{} _ = return (coerce x) -power _ y@QNaN{} = return (coerce y) - -powerSign :: Decimal a b -> Decimal c d -> Sign -powerSign x y - | Number.isNegative x && fromMaybe False (odd <$> integralValue y) = Neg - | otherwise = Pos - -integralPower :: (Precision p, Rounding r) - => Decimal a b -> Integer -> Arith p r (Decimal p r) -integralPower b e = integralPower' (return b) e one - where integralPower' :: (Precision p, Rounding r) - => Arith p r (Decimal a b) -> Integer -> Decimal p r - -> Arith p r (Decimal p r) - integralPower' _ 0 r = return r - integralPower' mb e r - | odd e = mb >>= \b -> multiply r b >>= - integralPower' (multiply b b) e' - | otherwise = integralPower' (mb >>= \b -> multiply b b) e' r - where e' = e `div` 2 - -{- $doctest-power ->>> op2 Op.power "2" "3" -8 - ->>> op2 Op.power "-2" "3" --8 - ->>> op2 Op.power "2" "-3" -0.125 - ->>> op2 Op.power "1.7" "8" -69.7575744 - ->>> op2 Op.power "10" "0.301029996" -2.00000000 - ->>> op2 Op.power "Infinity" "-1" -0 - ->>> op2 Op.power "Infinity" "0" -1 - ->>> op2 Op.power "Infinity" "1" -Infinity - ->>> op2 Op.power "-Infinity" "-1" --0 - ->>> op2 Op.power "-Infinity" "0" -1 - ->>> op2 Op.power "-Infinity" "1" --Infinity - ->>> op2 Op.power "-Infinity" "2" -Infinity - ->>> op2 Op.power "0" "0" -NaN --} - --- | 'quantize' takes two operands. If either operand is a /special value/ --- then the general rules apply, except that if either operand is infinite and --- the other is finite an Invalid operation condition is raised and the result --- is NaN, or if both are infinite then the result is the first operand. --- --- Otherwise (both operands are finite), 'quantize' returns the number which --- is equal in value (except for any rounding) and sign to the first --- (left-hand) operand and which has an /exponent/ set to be equal to the --- exponent of the second (right-hand) operand. --- --- The /coefficient/ of the result is derived from that of the left-hand --- operand. It may be rounded using the current /rounding/ setting (if the --- /exponent/ is being increased), multiplied by a positive power of ten (if --- the /exponent/ is being decreased), or is unchanged (if the /exponent/ is --- already equal to that of the right-hand operand). --- --- Unlike other operations, if the length of the /coefficient/ after the --- quantize operation would be greater than /precision/ then an Invalid --- operation condition is raised. This guarantees that, unless there is an --- error condition, the /exponent/ of the result of a quantize is always equal --- to that of the right-hand operand. --- --- Also unlike other operations, quantize will never raise Underflow, even if --- the result is subnormal and inexact. -quantize :: (Precision p, Rounding r) - => Decimal p r -> Decimal a b -> Arith p r (Decimal p r) -quantize x@Num { coefficient = xc, exponent = xe } Num { exponent = ye } - | xe > ye = result x { coefficient = xc * 10^(xe - ye), exponent = ye } - | xe < ye = rc >>= \c -> return x { coefficient = c, exponent = ye } - | otherwise = return x - - where result :: Precision p => Decimal p r -> Arith p r (Decimal p r) - result x = getPrecision >>= \p -> case numDigits (coefficient x) of - n | maybe False (n >) p -> invalidOperation x - _ -> return x - - rc :: Rounding r => Arith p r Coefficient - rc = let b = 10^(ye - xe) - (q, r) = xc `quotRem` b - in getRounder >>= \rounder -> return (rounder (sign x) r b q) - -quantize Num{} Inf{} = invalidOperation qNaN -quantize Inf{} Num{} = invalidOperation qNaN -quantize n@Inf{} Inf{} = return n -quantize n@SNaN{} _ = invalidOperation n -quantize _ n@SNaN{} = invalidOperation n -quantize n@QNaN{} _ = return n -quantize _ n@QNaN{} = return (coerce n) - -{- $doctest-quantize ->>> op2 Op.quantize "2.17" "0.001" -2.170 - ->>> op2 Op.quantize "2.17" "0.01" -2.17 - ->>> op2 Op.quantize "2.17" "0.1" -2.2 - ->>> op2 Op.quantize "2.17" "1e+0" -2 - ->>> op2 Op.quantize "2.17" "1e+1" -0E+1 - ->>> op2 Op.quantize "-Inf" "Infinity" --Infinity - ->>> op2 Op.quantize "2" "Infinity" -NaN - ->>> op2 Op.quantize "-0.1" "1" --0 - ->>> op2 Op.quantize "-0" "1e+5" --0E+5 - ->>> op2 Op.quantize "+35236450.6" "1e-2" -NaN - ->>> op2 Op.quantize "-35236450.6" "1e-2" -NaN - ->>> op2 Op.quantize "217" "1e-1" -217.0 - ->>> op2 Op.quantize "217" "1e+0" -217 - ->>> op2 Op.quantize "217" "1e+1" -2.2E+2 - ->>> op2 Op.quantize "217" "1e+2" -2E+2 --} - --- | 'reduce' takes one operand. It has the same semantics as the 'plus' --- operation, except that if the final result is finite it is reduced to its --- simplest form, with all trailing zeros removed and its sign preserved. -reduce :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r) -reduce n = reduce' <$> plus n - where reduce' n@Num { coefficient = c, exponent = e } - | c == 0 = n { exponent = 0 } - | r == 0 = reduce' n { coefficient = q, exponent = e + 1 } - where (q, r) = c `quotRem` 10 - reduce' n = n - -{- $doctest-reduce ->>> op1 Op.reduce "2.1" -2.1 - ->>> op1 Op.reduce "-2.0" --2 - ->>> op1 Op.reduce "1.200" -1.2 - ->>> op1 Op.reduce "-120" --1.2E+2 - ->>> op1 Op.reduce "120.00" -1.2E+2 - ->>> op1 Op.reduce "0.00" -0 --} - --- | 'squareRoot' takes one operand. If the operand is a /special value/ then --- the general rules apply. --- --- Otherwise, the ideal exponent of the result is defined to be half the --- exponent of the operand (rounded to an integer, towards −Infinity, if --- necessary) and then: --- --- If the operand is less than zero an Invalid operation condition is raised. --- --- If the operand is greater than zero, the result is the square root of the --- operand. If no rounding is necessary (the exact result requires /precision/ --- digits or fewer) then the coefficient and exponent giving the correct value --- and with the exponent closest to the ideal exponent is used. If the result --- must be inexact, it is rounded using the /round-half-even/ algorithm and --- the coefficient will have exactly /precision/ digits (unless the result is --- subnormal), and the exponent will be set to maintain the correct value. --- --- Otherwise (the operand is equal to zero), the result will be the zero with --- the same sign as the operand and with the ideal exponent. -squareRoot :: FinitePrecision p - => Decimal a b -> Arith p r (Decimal p RoundHalfEven) -squareRoot n@Num { sign = s, coefficient = c, exponent = e } - | c == 0 = return n { exponent = idealExp } - | s == Pos = subResult >>= subRounded >>= result - - where idealExp = e `div` 2 :: Exponent - - reduced :: Decimal p r -> Decimal p r - reduced n@Num { coefficient = c, exponent = e } - | e < idealExp = case bd of - Just (b, (q, _)) -> n { coefficient = q, exponent = e + b } - Nothing -> n - | e > idealExp = n { coefficient = c * 10^d, exponent = idealExp } - where d = Prelude.abs (e - idealExp) - bd = find (\(_, (_, r)) -> r == 0) ds - ds = map (\d -> (d, c `quotRem` (10^d))) [d, d - 1 .. 1] - reduced n = n - - subResult :: FinitePrecision p - => Arith p r (Decimal (PPlus1 (PPlus1 p)) RoundHalfEven) - subResult = subArith (babylonian approx) - - subRounded :: Precision p - => Decimal a b -> Arith p r (Decimal p RoundHalfEven) - subRounded = subArith . roundDecimal - - exactness :: Decimal a b -> Arith p r (Decimal PInfinite RoundHalfEven) - exactness r = subArith (multiply r r >>= compare n) - - result :: Decimal p a -> Arith p r (Decimal p a) - result r = do - e <- exactness r - if Number.isZero e - then return (reduced r) - else let r' = coerce r - in coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r') - - approx :: Decimal p r - approx | even ae = n { coefficient = 2, exponent = ae `quot` 2 } - | otherwise = n { coefficient = 6, exponent = (ae - 1) `quot` 2 } - where ae = adjustedExponent n - - babylonian :: FinitePrecision p => Decimal p RoundHalfEven - -> Arith p RoundHalfEven (Decimal p RoundHalfEven) - babylonian x = do - x' <- multiply oneHalf =<< add x =<< n `divide` x - if x' == x then return x' else babylonian x' - -squareRoot n@Inf { sign = Pos } = return (coerce n) -squareRoot n@QNaN{} = return (coerce n) -squareRoot n = coerce <$> invalidOperation n - -{- $doctest-squareRoot ->>> op1 Op.squareRoot "0" -0 - ->>> op1 Op.squareRoot "-0" --0 - -This example appears to contradict the specification that the resulting -coefficient will have exactly /precision/ digits; awaiting clarification. -<<< op1 Op.squareRoot "0.39" -0.62449980 - ->>> op1 Op.squareRoot "100" -10 - ->>> op1 Op.squareRoot "1" -1 - ->>> op1 Op.squareRoot "1.0" -1.0 - ->>> op1 Op.squareRoot "1.00" -1.0 - ->>> op1 Op.squareRoot "7" -2.64575131 - ->>> op1 Op.squareRoot "10" -3.16227766 --} - --- $miscellaneous-operations --- --- This section describes miscellaneous operations on decimal numbers, --- including non-numeric comparisons, sign and other manipulations, and --- logical operations. --- --- The logical operations ('and', 'invert', 'or', and 'xor') take --- /logical operands/, which are finite numbers with a /sign/ of 0, an --- /exponent/ of 0, and a /coefficient/ whose digits must all be either 0 or --- 1. The length of the result will be at most /precision/ digits (all of --- which will be either 0 or 1); operands are truncated on the left or padded --- with zeros on the left as necessary. The result of a logical operation is --- never rounded and the only /flag/ that might be set is /invalid-operation/ --- (set if an operand is not a valid logical operand). --- --- Some operations return a boolean value that is described as 0 or 1 in the --- documentation below. For reasons of efficiency, and as permitted by the --- /General Decimal Arithmetic Specification/, these operations return a --- 'Bool' in this implementation, but can be converted to 'Decimal' via --- 'fromBool'. - -data Logical = Logical { bits :: Integer, bitLength :: Int } - -toLogical :: Decimal a b -> Maybe Logical -toLogical Num { sign = Pos, coefficient = c, exponent = 0 } = - getBits c Logical { bits = zeroBits, bitLength = 0 } - - where getBits :: Coefficient -> Logical -> Maybe Logical - getBits 0 g = return g - getBits c g@Logical { bits = b, bitLength = l } = case d of - 0 -> getBits c' g { bitLength = succ l } - 1 -> getBits c' g { bits = setBit b l, bitLength = succ l } - _ -> Nothing - where (c', d) = c `quotRem` 10 - -toLogical _ = Nothing - -fromLogical :: Logical -> Decimal a b -fromLogical Logical { bits = b, bitLength = l } = - Num { sign = Pos, coefficient = fromBits 0 1 0, exponent = 0 } - - where fromBits :: Int -> Coefficient -> Coefficient -> Coefficient - fromBits i r c - | i == l = c - | testBit b i = fromBits i' r' (c + r) - | otherwise = fromBits i' r' c - where i' = succ i - r' = r * 10 - --- | 'and' is a logical operation which takes two logical operands. The result --- is the digit-wise /and/ of the two operands; each digit of the result is --- the logical and of the corresponding digits of the operands, aligned at the --- least-significant digit. A result digit is 1 if both of the corresponding --- operand digits are 1; otherwise it is 0. -and :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) -and x y = case (toLogical x, toLogical y) of - (Just lx, Just ly) -> getPrecision >>= \p -> - let m = Prelude.min (bitLength lx) (bitLength ly) - z = Logical { bits = bits lx .&. bits ly - , bitLength = maybe m (Prelude.min m) p } - in return (fromLogical z) - _ -> invalidOperation qNaN - -{- $doctest-and ->>> op2 Op.and "0" "0" -0 - ->>> op2 Op.and "0" "1" -0 - ->>> op2 Op.and "1" "0" -0 - ->>> op2 Op.and "1" "1" -1 - ->>> op2 Op.and "1100" "1010" -1000 - ->>> op2 Op.and "1111" "10" -10 --} - --- | 'or' is a logical operation which takes two logical operands. The result --- is the digit-wise /inclusive or/ of the two operands; each digit of the --- result is the logical or of the corresponding digits of the operands, --- aligned at the least-significant digit. A result digit is 1 if either or --- both of the corresponding operand digits is 1; otherwise it is 0. -or :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) -or x y = case (toLogical x, toLogical y) of - (Just lx, Just ly) -> getPrecision >>= \p -> - let m = Prelude.max (bitLength lx) (bitLength ly) - z = Logical { bits = bits lx .|. bits ly - , bitLength = maybe m (Prelude.min m) p } - in return (fromLogical z) - _ -> invalidOperation qNaN - -{- $doctest-or ->>> op2 Op.or "0" "0" -0 - ->>> op2 Op.or "0" "1" -1 - ->>> op2 Op.or "1" "0" -1 - ->>> op2 Op.or "1" "1" -1 - ->>> op2 Op.or "1100" "1010" -1110 - ->>> op2 Op.or "1110" "10" -1110 --} - --- | 'xor' is a logical operation which takes two logical operands. The result --- is the digit-wise /exclusive or/ of the two operands; each digit of the --- result is the logical exclusive-or of the corresponding digits of the --- operands, aligned at the least-significant digit. A result digit is 1 if --- one of the corresponding operand digits is 1 and the other is 0; otherwise --- it is 0. -xor :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) -xor x y = case (toLogical x, toLogical y) of - (Just lx, Just ly) -> getPrecision >>= \p -> - let m = Prelude.max (bitLength lx) (bitLength ly) - z = Logical { bits = bits lx `Bits.xor` bits ly - , bitLength = maybe m (Prelude.min m) p } - in return (fromLogical z) - _ -> invalidOperation qNaN - -{- $doctest-xor ->>> op2 Op.xor "0" "0" -0 - ->>> op2 Op.xor "0" "1" -1 - ->>> op2 Op.xor "1" "0" -1 - ->>> op2 Op.xor "1" "1" -0 - ->>> op2 Op.xor "1100" "1010" -110 - ->>> op2 Op.xor "1111" "10" -1101 --} - --- | 'invert' is a logical operation which takes one logical operand. The --- result is the digit-wise /inversion/ of the operand; each digit of the --- result is the inverse of the corresponding digit of the operand. A result --- digit is 1 if the corresponding operand digit is 0; otherwise it is 0. -invert :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p r) -invert x = case toLogical x of - Just lx -> getPrecision >>= \(Just p) -> - let z = Logical { bits = complement (bits lx), bitLength = p } - in return (fromLogical z) - _ -> invalidOperation qNaN - -{- $doctest-invert ->>> op1 Op.invert "0" -111111111 - ->>> op1 Op.invert "1" -111111110 - ->>> op1 Op.invert "111111111" -0 - ->>> op1 Op.invert "101010101" -10101010 --} - --- | 'canonical' takes one operand. The result has the same value as the --- operand but always uses a /canonical/ encoding. The definition of --- /canonical/ is implementation-defined; if more than one internal encoding --- for a given NaN, Infinity, or finite number is possible then one --- “preferred” encoding is deemed canonical. This operation then returns the --- value using that preferred encoding. --- --- If all possible operands have just one internal encoding each, then --- 'canonical' always returns the operand unchanged (that is, it has the same --- effect as 'copy'). This operation is unaffected by context and is quiet — --- no /flags/ are changed in the context. -canonical :: Decimal a b -> Arith p r (Decimal a b) -canonical = return - -{- $doctest-canonical ->>> op1 Op.canonical "2.50" -2.50 --} - --- | 'class_' takes one operand. The result is an indication of the /class/ of --- the operand, where the class is one of ten possibilities, corresponding to --- one of the strings @"sNaN"@ (signaling NaN), @\"NaN"@ (quiet NaN), --- @"-Infinity"@ (negative infinity), @"-Normal"@ (negative normal finite --- number), @"-Subnormal"@ (negative subnormal finite number), @"-Zero"@ --- (negative zero), @"+Zero"@ (non-negative zero), @"+Subnormal"@ (positive --- subnormal finite number), @"+Normal"@ (positive normal finite number), or --- @"+Infinity"@ (positive infinity). This operation is quiet; no /flags/ are --- changed in the context. --- --- Note that unlike the special values in the model, the sign of any NaN is --- ignored in the classification, as required by IEEE 754. -class_ :: Precision a => Decimal a b -> Arith p r Class -class_ n = return $ case n of - Num {} | Number.isZero n -> NumberClass (sign n) ZeroClass - | Number.isSubnormal n -> NumberClass (sign n) SubnormalClass - | otherwise -> NumberClass (sign n) NormalClass - Inf {} -> NumberClass (sign n) InfinityClass - QNaN{} -> NaNClass QNaNClass - SNaN{} -> NaNClass SNaNClass - -data Class = NumberClass Sign NumberClass -- ^ Number (finite or infinite) - | NaNClass NaNClass -- ^ Not a number (quiet or signaling) - deriving Eq - -data NumberClass = ZeroClass -- ^ Zero - | SubnormalClass -- ^ Subnormal finite number - | NormalClass -- ^ Normal finite number - | InfinityClass -- ^ Infinity - deriving Eq - -data NaNClass = QNaNClass -- ^ Not a number (quiet) - | SNaNClass -- ^ Not a number (signaling) - deriving Eq - -instance Show Class where - show c = case c of - NumberClass s nc -> signChar s : showNumberClass nc - NaNClass QNaNClass -> nan - NaNClass SNaNClass -> 's' : nan - - where signChar :: Sign -> Char - signChar Pos = '+' - signChar Neg = '-' - - showNumberClass :: NumberClass -> String - showNumberClass s = case s of - ZeroClass -> "Zero" - SubnormalClass -> "Subnormal" - NormalClass -> "Normal" - InfinityClass -> "Infinity" - - nan :: String - nan = "NaN" - -{- $doctest-class_ ->>> op1 Op.class_ "Infinity" -+Infinity - ->>> op1 Op.class_ "1E-10" -+Normal - ->>> op1 Op.class_ "2.50" -+Normal - ->>> op1 Op.class_ "0.1E-999" -+Subnormal - ->>> op1 Op.class_ "0" -+Zero - ->>> op1 Op.class_ "-0" --Zero - ->>> op1 Op.class_ "-0.1E-999" --Subnormal - ->>> op1 Op.class_ "-1E-10" --Normal - ->>> op1 Op.class_ "-2.50" --Normal - ->>> op1 Op.class_ "-Infinity" --Infinity - ->>> op1 Op.class_ "NaN" -NaN - ->>> op1 Op.class_ "-NaN" -NaN - ->>> op1 Op.class_ "sNaN" -sNaN --} - --- | 'copy' takes one operand. The result is a copy of the operand. This --- operation is unaffected by context and is quiet — no /flags/ are changed in --- the context. -copy :: Decimal a b -> Arith p r (Decimal a b) -copy = return - -{- $doctest-copy ->>> op1 Op.copy "2.1" -2.1 - ->>> op1 Op.copy "-1.00" --1.00 --} - --- | 'copyAbs' takes one operand. The result is a copy of the operand with the --- /sign/ set to 0. Unlike the 'abs' operation, this operation is unaffected --- by context and is quiet — no /flags/ are changed in the context. -copyAbs :: Decimal a b -> Arith p r (Decimal a b) -copyAbs n = return n { sign = Pos } - -{- $doctest-copyAbs ->>> op1 Op.copyAbs "2.1" -2.1 - ->>> op1 Op.copyAbs "-100" -100 --} - --- | 'copyNegate' takes one operand. The result is a copy of the operand with --- the /sign/ inverted (a /sign/ of 0 becomes 1 and vice versa). Unlike the --- 'minus' operation, this operation is unaffected by context and is quiet — --- no /flags/ are changed in the context. -copyNegate :: Decimal a b -> Arith p r (Decimal a b) -copyNegate n = return n { sign = negateSign (sign n) } - -{- $doctest-copyNegate ->>> op1 Op.copyNegate "101.5" --101.5 - ->>> op1 Op.copyNegate "-101.5" -101.5 --} - --- | 'copySign' takes two operands. The result is a copy of the first operand --- with the /sign/ set to be the same as the /sign/ of the second --- operand. This operation is unaffected by context and is quiet — no /flags/ --- are changed in the context. -copySign :: Decimal a b -> Decimal c d -> Arith p r (Decimal a b) -copySign n m = return n { sign = sign m } - -{- $doctest-copySign ->>> op2 Op.copySign "1.50" "7.33" -1.50 - ->>> op2 Op.copySign "-1.50" "7.33" -1.50 - ->>> op2 Op.copySign "1.50" "-7.33" --1.50 - ->>> op2 Op.copySign "-1.50" "-7.33" --1.50 --} - --- | 'isCanonical' takes one operand. The result is 1 if the operand is --- /canonical/; otherwise it is 0. The definition of /canonical/ is --- implementation-defined; if more than one internal encoding for a given NaN, --- Infinity, or finite number is possible then one “preferred” encoding is --- deemed canonical. This operation then tests whether the internal encoding --- is that preferred encoding. --- --- If all possible operands have just one internal encoding each, then --- 'isCanonical' always returns 1. This operation is unaffected by context and --- is quiet — no /flags/ are changed in the context. -isCanonical :: Decimal a b -> Arith p r Bool -isCanonical _ = return True - -{- $doctest-isCanonical ->>> fromBool $ op1 Op.isCanonical "2.50" -1 --} - --- | 'isFinite' takes one operand. The result is 1 if the operand is neither --- infinite nor a NaN (that is, it is a normal number, a subnormal number, or --- a zero); otherwise it is 0. This operation is unaffected by context and is --- quiet — no /flags/ are changed in the context. -isFinite :: Decimal a b -> Arith p r Bool -isFinite = return . Number.isFinite - -{- $doctest-isFinite ->>> fromBool $ op1 Op.isFinite "2.50" -1 - ->>> fromBool $ op1 Op.isFinite "-0.3" -1 - ->>> fromBool $ op1 Op.isFinite "0" -1 - ->>> fromBool $ op1 Op.isFinite "Inf" -0 - ->>> fromBool $ op1 Op.isFinite "NaN" -0 --} - --- | 'isInfinite' takes one operand. The result is 1 if the operand is an --- Infinity; otherwise it is 0. This operation is unaffected by context and is --- quiet — no /flags/ are changed in the context. -isInfinite :: Decimal a b -> Arith p r Bool -isInfinite n = return $ case n of - Inf{} -> True - _ -> False - -{- $doctest-isInfinite ->>> fromBool $ op1 Op.isInfinite "2.50" -0 - ->>> fromBool $ op1 Op.isInfinite "-Inf" -1 - ->>> fromBool $ op1 Op.isInfinite "NaN" -0 --} - --- | 'isNaN' takes one operand. The result is 1 if the operand is a NaN (quiet --- or signaling); otherwise it is 0. This operation is unaffected by context --- and is quiet — no /flags/ are changed in the context. -isNaN :: Decimal a b -> Arith p r Bool -isNaN n = return $ case n of - QNaN{} -> True - SNaN{} -> True - _ -> False - -{- $doctest-isNaN ->>> fromBool $ op1 Op.isNaN "2.50" -0 - ->>> fromBool $ op1 Op.isNaN "NaN" -1 - ->>> fromBool $ op1 Op.isNaN "-sNaN" -1 --} - --- | 'isNormal' takes one operand. The result is 1 if the operand is a --- positive or negative /normal number/; otherwise it is 0. This operation is --- quiet; no /flags/ are changed in the context. -isNormal :: Precision a => Decimal a b -> Arith p r Bool -isNormal = return . Number.isNormal - -{- $doctest-isNormal ->>> fromBool $ op1 Op.isNormal "2.50" -1 - ->>> fromBool $ op1 Op.isNormal "0.1E-999" -0 - ->>> fromBool $ op1 Op.isNormal "0.00" -0 - ->>> fromBool $ op1 Op.isNormal "-Inf" -0 - ->>> fromBool $ op1 Op.isNormal "NaN" -0 --} - --- | 'isQNaN' takes one operand. The result is 1 if the operand is a quiet --- NaN; otherwise it is 0. This operation is unaffected by context and is --- quiet — no /flags/ are changed in the context. -isQNaN :: Decimal a b -> Arith p r Bool -isQNaN n = return $ case n of - QNaN{} -> True - _ -> False - -{- $doctest-isQNaN ->>> fromBool $ op1 Op.isQNaN "2.50" -0 - ->>> fromBool $ op1 Op.isQNaN "NaN" -1 - ->>> fromBool $ op1 Op.isQNaN "sNaN" -0 --} - --- | 'isSigned' takes one operand. The result is 1 if the /sign/ of the --- operand is 1; otherwise it is 0. This operation is unaffected by context --- and is quiet — no /flags/ are changed in the context. -isSigned :: Decimal a b -> Arith p r Bool -isSigned = return . Number.isNegative - -{- $doctest-isSigned ->>> fromBool $ op1 Op.isSigned "2.50" -0 - ->>> fromBool $ op1 Op.isSigned "-12" -1 - ->>> fromBool $ op1 Op.isSigned "-0" -1 --} - --- | 'isSNaN' takes one operand. The result is 1 if the operand is a signaling --- NaN; otherwise it is 0. This operation is unaffected by context and is --- quiet — no /flags/ are changed in the context. -isSNaN :: Decimal a b -> Arith p r Bool -isSNaN n = return $ case n of - SNaN{} -> True - _ -> False - -{- $doctest-isSNaN ->>> fromBool $ op1 Op.isSNaN "2.50" -0 - ->>> fromBool $ op1 Op.isSNaN "NaN" -0 - ->>> fromBool $ op1 Op.isSNaN "sNaN" -1 --} - --- | 'isSubnormal' takes one operand. The result is 1 if the operand is a --- positive or negative /subnormal number/; otherwise it is 0. This operation --- is quiet; no /flags/ are changed in the context. -isSubnormal :: Precision a => Decimal a b -> Arith p r Bool -isSubnormal = return . Number.isSubnormal - -{- $doctest-isSubnormal ->>> fromBool $ op1 Op.isSubnormal "2.50" -0 - ->>> fromBool $ op1 Op.isSubnormal "0.1E-999" -1 - ->>> fromBool $ op1 Op.isSubnormal "0.00" -0 - ->>> fromBool $ op1 Op.isSubnormal "-Inf" -0 - ->>> fromBool $ op1 Op.isSubnormal "NaN" -0 --} - --- | 'isZero' takes one operand. The result is 1 if the operand is a zero; --- otherwise it is 0. This operation is unaffected by context and is quiet — --- no /flags/ are changed in the context. -isZero :: Decimal a b -> Arith p r Bool -isZero = return . Number.isZero - -{- $doctest-isZero ->>> fromBool $ op1 Op.isZero "0" -1 - ->>> fromBool $ op1 Op.isZero "2.50" -0 - ->>> fromBool $ op1 Op.isZero "-0E+2" -1 --} - --- | 'logb' takes one operand. If the operand is a NaN then the general --- arithmetic rules apply. If the operand is infinite then +Infinity is --- returned. If the operand is a zero, then −Infinity is returned and the --- Division by zero exceptional condition is raised. --- --- Otherwise, the result is the integer which is the exponent of the magnitude --- of the most significant digit of the operand (as though the operand were --- truncated to a single digit while maintaining the value of that digit and --- without limiting the resulting exponent). All results are exact unless an --- integer result does not fit in the available /precision/. -logb :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r) -logb Num { coefficient = c, exponent = e } - | c == 0 = raiseSignal DivisionByZero Inf { sign = Neg } - | otherwise = roundDecimal (fromInteger r :: Decimal PInfinite RoundHalfEven) - where r = fromIntegral (numDigits c) - 1 + fromIntegral e :: Integer -logb Inf{} = return Inf { sign = Pos } -logb n@QNaN{} = return (coerce n) -logb n@SNaN{} = invalidOperation n - -{- $doctest-logb ->>> op1 Op.logb "250" -2 - ->>> op1 Op.logb "2.50" -0 - ->>> op1 Op.logb "0.03" --2 - ->>> op1 Op.logb "0" --Infinity --} - --- | 'radix' takes no operands. The result is the radix (base) in which --- arithmetic is effected; for this specification the result will have the --- value 10. -radix :: Precision p => Arith p r (Decimal p r) -radix = return radix' - where radix' = case precision radix' of - Just 1 -> one { exponent = 1 } - _ -> one { coefficient = 10 } - -{- $doctest-radix ->>> op0 Op.radix -10 --} - --- | 'sameQuantum' takes two operands, and returns 1 if the two operands have --- the same /exponent/ or 0 otherwise. The result is never affected by either --- the sign or the coefficient of either operand. --- --- If either operand is a /special value/, 1 is returned only if both operands --- are NaNs or both are infinities. --- --- 'sameQuantum' does not change any /flags/ in the context. -sameQuantum :: Decimal a b -> Decimal c d -> Arith p r Bool -sameQuantum Num { exponent = e1 } Num { exponent = e2 } - | e1 == e2 = return True - | otherwise = return False -sameQuantum Inf {} Inf {} = return True -sameQuantum QNaN{} QNaN{} = return True -sameQuantum SNaN{} SNaN{} = return True -sameQuantum QNaN{} SNaN{} = return True -sameQuantum SNaN{} QNaN{} = return True -sameQuantum _ _ = return False - -{- $doctest-sameQuantum ->>> fromBool $ op2 Op.sameQuantum "2.17" "0.001" -0 - ->>> fromBool $ op2 Op.sameQuantum "2.17" "0.01" -1 - ->>> fromBool $ op2 Op.sameQuantum "2.17" "0.1" -0 - ->>> fromBool $ op2 Op.sameQuantum "2.17" "1" -0 - ->>> fromBool $ op2 Op.sameQuantum "Inf" "-Inf" -1 - ->>> fromBool $ op2 Op.sameQuantum "NaN" "NaN" -1 --} - --- | 'shift' takes two operands. The second operand must be an integer (with --- an /exponent/ of 0) in the range /−precision/ through /precision/. If the --- first operand is a NaN then the general arithmetic rules apply, and if it --- is infinite then the result is the Infinity unchanged. --- --- Otherwise (the first operand is finite) the result has the same /sign/ and --- /exponent/ as the first operand, and a /coefficient/ which is a shifted --- copy of the digits in the coefficient of the first operand. The number of --- places to shift is taken from the absolute value of the second operand, --- with the shift being to the left if the second operand is positive or to --- the right otherwise. Digits shifted into the coefficient are zeros. --- --- The only /flag/ that might be set is /invalid-operation/ (set if the first --- operand is an sNaN or the second is not valid). --- --- The 'rotate' operation can be used to rotate rather than shift a --- coefficient. -shift :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) -shift n@Num { coefficient = c } s@Num { sign = d, coefficient = sc } - | validShift z s = return z - where z = case precision z of - Just p -> y { coefficient = coefficient y `rem` 10 ^ p } - Nothing -> y - y = case d of - Pos -> n { coefficient = c * 10 ^ sc } - Neg -> n { coefficient = c `quot` 10 ^ sc } -shift n@Inf{} s | validShift z s = return z where z = coerce n -shift n@QNaN{} s | validShift z s = return z where z = coerce n -shift n _ = invalidOperation n - -validShift :: Precision p => p -> Decimal a b -> Bool -validShift px Num { coefficient = c, exponent = 0 } = - let p = fromIntegral <$> precision px in maybe True (c <=) p -validShift _ _ = False - -{- $doctest-shift ->>> op2 Op.shift "34" "8" -400000000 - ->>> op2 Op.shift "12" "9" -0 - ->>> op2 Op.shift "123456789" "-2" -1234567 - ->>> op2 Op.shift "123456789" "0" -123456789 - ->>> op2 Op.shift "123456789" "+2" -345678900 --} - --- | 'rotate' takes two operands. The second operand must be an integer (with --- an /exponent/ of 0) in the range /−precision/ through /precision/. If the --- first operand is a NaN then the general arithmetic rules apply, and if it --- is infinite then the result is the Infinity unchanged. --- --- Otherwise (the first operand is finite) the result has the same /sign/ and --- /exponent/ as the first operand, and a /coefficient/ which is a rotated --- copy of the digits in the coefficient of the first operand. The number of --- places of rotation is taken from the absolute value of the second operand, --- with the rotation being to the left if the second operand is positive or to --- the right otherwise. --- --- If the coefficient of the first operand has fewer than /precision/ digits, --- it is treated as though it were padded on the left with zeros to length --- /precision/ before the rotation. Similarly, if the coefficient of the first --- operand has more than /precision/ digits, it is truncated on the left --- before use. --- --- The only /flag/ that might be set is /invalid-operation/ (set if the first --- operand is an sNaN or the second is not valid). --- --- The 'shift' operation can be used to shift rather than rotate a --- coefficient. -rotate :: FinitePrecision p - => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) -rotate n@Num { coefficient = c } s@Num { sign = d, coefficient = sc } - | validShift z s = return z - where z = n { coefficient = rc * b + (lc `rem` b) } - (lc, rc) = c `quotRem` b' - (b , b') = case d of - Pos -> (10^sc , 10^sc') - Neg -> (10^sc', 10^sc ) - Just p = precision z - sc' = p - fromIntegral sc -rotate n@Inf{} s | validShift z s = return z where z = coerce n -rotate n@QNaN{} s | validShift z s = return z where z = coerce n -rotate n _ = invalidOperation n - -{- $doctest-rotate ->>> op2 Op.rotate "34" "8" -400000003 - ->>> op2 Op.rotate "12" "9" -12 - ->>> op2 Op.rotate "123456789" "-2" -891234567 - ->>> op2 Op.rotate "123456789" "0" -123456789 - ->>> op2 Op.rotate "123456789" "+2" -345678912 --} + -- ** General arithmetic + add + , subtract + , multiply + , divide + -- divideInteger + -- remainder + -- remainderNear + , power + , squareRoot + , fusedMultiplyAdd + + -- ** Exponential and logarithmic + , exp + , ln + , log10 + + -- ** Unary sign + , plus + , minus + , abs + + -- ** Comparison + , compare + , compareSignal + + , min + , max + , minMagnitude + , maxMagnitude + + -- ** Rounding and quantization + + , roundToIntegralValue + , roundToIntegralExact + , quantize + , reduce + + -- nextMinus + -- nextPlus + -- nextToward + + -- * Miscellaneous operations + -- $miscellaneous-operations + + -- ** Logic and shifting + -- $logical-operations + , and + , or + , xor + , invert + + , shift + , rotate + + -- ** Predicates + , isZero + , isSigned + , isFinite + , isInfinite + , isNormal + , isSubnormal + , isNaN + , isQNaN + , isSNaN + , isCanonical + + -- ** Total comparison and classification + , compareTotal + , compareTotalMagnitude + + , class', Class(..), Sign(..), NumberClass(..), NaNClass(..) + + -- ** Exponent manipulation + , logb + , scaleb + , sameQuantum + , radix + + -- ** Sign manipulation and conversion + , copyAbs + , copyNegate + , copySign + , copy + + , canonical + ) where + +import Prelude hiding (abs, and, compare, exp, exponent, isInfinite, isNaN, + max, min, or, subtract) +import qualified Prelude + +import Control.Monad (join) +import Data.Bits (complement, setBit, testBit, zeroBits, (.&.), (.|.)) +import Data.Coerce (coerce) +import Data.List (find) +import Data.Maybe (fromMaybe) + +import qualified Data.Bits as Bits + +import Numeric.Decimal.Arithmetic +import Numeric.Decimal.Exception +import Numeric.Decimal.Number hiding (isFinite, isNormal, isSubnormal, isZero) +import Numeric.Decimal.Precision +import Numeric.Decimal.Rounding + +import qualified Numeric.Decimal.Number as Number + +finitePrecision :: FinitePrecision p => Decimal p r -> Int +finitePrecision n = let Just p = precision n in p + +roundingAlg :: Rounding r => Arith p r a -> RoundingAlgorithm +roundingAlg = rounding . arithRounding + where arithRounding :: Arith p r a -> r + arithRounding = undefined + +result :: (Precision p, Rounding r) => Decimal p r -> Arith p r (Decimal p r) +result = roundDecimal -- ... +-- | maybe False (numDigits c >) (precision r) = undefined + +generalRules1 :: Decimal a b -> Arith p r (Decimal p r) +generalRules1 nan@NaN { signaling = False } = return (coerce nan) +generalRules1 x = invalidOperation x + +generalRules2 :: Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +generalRules2 nan@NaN { signaling = True } _ = invalidOperation nan +generalRules2 _ nan@NaN { signaling = True } = invalidOperation nan +generalRules2 nan@NaN{} _ = return (coerce nan) +generalRules2 _ nan@NaN{} = return (coerce nan) +generalRules2 x _ = invalidOperation x + +-- $arithmetic-operations +-- +-- This section describes the arithmetic operations on, and some other +-- functions of, numbers, including subnormal numbers, negative zeros, and +-- special values (see also IEEE 754 §5 and §6). + +-- | 'add' takes two operands. If either operand is a /special value/ then the +-- general rules apply. +-- +-- Otherwise, the operands are added. +-- +-- The result is then rounded to /precision/ digits if necessary, counting +-- from the most significant digit of the result. +add :: (Precision p, Rounding r) + => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +add Num { sign = xs, coefficient = xc, exponent = xe } + Num { sign = ys, coefficient = yc, exponent = ye } = sum + + where sum = result Num { sign = rs, coefficient = rc, exponent = re } + rs | rc /= 0 = if xac > yac then xs else ys + | xs == Neg && ys == Neg = Neg + | xs /= ys && + roundingAlg sum == RoundFloor = Neg + | otherwise = Pos + rc | xs == ys = xac + yac + | xac > yac = xac - yac + | otherwise = yac - xac + re = Prelude.min xe ye + (xac, yac) | xe == ye = (xc, yc) + | xe > ye = (xc * 10^n, yc) + | otherwise = (xc, yc * 10^n) + where n = Prelude.abs (xe - ye) + +add inf@Inf{} Num{} = return (coerce inf) +add Num{} inf@Inf{} = return (coerce inf) +add inf@Inf { sign = xs } Inf { sign = ys } + | xs == ys = return (coerce inf) + | otherwise = invalidOperation qNaN +add x y = generalRules2 x y + +-- | 'subtract' takes two operands. If either operand is a /special value/ +-- then the general rules apply. +-- +-- Otherwise, the operands are added after inverting the /sign/ used for the +-- second operand. +-- +-- The result is then rounded to /precision/ digits if necessary, counting +-- from the most significant digit of the result. +subtract :: (Precision p, Rounding r) + => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +subtract x@Num{} y@Num{} = add x (flipSign y) +subtract x@Inf{} y@Num{} = add x (flipSign y) +subtract x@Num{} y@Inf{} = add x (flipSign y) +subtract x@Inf{} y@Inf{} = add x (flipSign y) +subtract x y = generalRules2 x y + +-- | 'minus' takes one operand, and corresponds to the prefix minus operator +-- in programming languages. +-- +-- Note that the result of this operation is affected by context and may set +-- /flags/. The 'copyNegate' operation may be used instead of 'minus' if this +-- is not desired. +minus :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r) +minus x = zero { exponent = exponent x } `subtract` x + +-- | 'plus' takes one operand, and corresponds to the prefix plus operator in +-- programming languages. +-- +-- Note that the result of this operation is affected by context and may set +-- /flags/. +plus :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r) +plus x = zero { exponent = exponent x } `add` x + +-- | 'multiply' takes two operands. If either operand is a /special value/ +-- then the general rules apply. Otherwise, the operands are multiplied +-- together (“long multiplication”), resulting in a number which may be as +-- long as the sum of the lengths of the two operands. +-- +-- The result is then rounded to /precision/ digits if necessary, counting +-- from the most significant digit of the result. +multiply :: (Precision p, Rounding r) + => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +multiply Num { sign = xs, coefficient = xc, exponent = xe } + Num { sign = ys, coefficient = yc, exponent = ye } = result rn + + where rn = Num { sign = rs, coefficient = rc, exponent = re } + rs = xorSigns xs ys + rc = xc * yc + re = xe + ye + +multiply Inf { sign = xs } Inf { sign = ys } = + return Inf { sign = xorSigns xs ys } +multiply Inf { sign = xs } Num { sign = ys, coefficient = yc } + | yc == 0 = invalidOperation qNaN + | otherwise = return Inf { sign = xorSigns xs ys } +multiply Num { sign = xs, coefficient = xc } Inf { sign = ys } + | xc == 0 = invalidOperation qNaN + | otherwise = return Inf { sign = xorSigns xs ys } +multiply x y = generalRules2 x y + +-- | 'exp' takes one operand. If the operand is a NaN then the general rules +-- for special values apply. +-- +-- Otherwise, the result is /e/ raised to the power of the operand, with the +-- following cases: +-- +-- * If the operand is −Infinity, the result is 0 and exact. +-- +-- * If the operand is a zero, the result is 1 and exact. +-- +-- * If the operand is +Infinity, the result is +Infinity and exact. +-- +-- * Otherwise the result is inexact and will be rounded using the +-- 'RoundHalfEven' algorithm. The coefficient will have exactly /precision/ +-- digits (unless the result is subnormal). These inexact results should be +-- correctly rounded, but may be up to 1 ulp (unit in last place) in error. +exp :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p RoundHalfEven) +exp x@Num { sign = s, coefficient = c } + | c == 0 = return one + | s == Neg = subArith (maclaurin x { sign = Pos } >>= reciprocal) >>= + subRounded >>= result + | otherwise = subArith (maclaurin x) >>= subRounded >>= result + + where multiplyExact :: Decimal a b -> Decimal c d + -> Arith PInfinite RoundHalfEven + (Decimal PInfinite RoundHalfEven) + multiplyExact = multiply + + maclaurin :: FinitePrecision p => Decimal a b + -> Arith p RoundHalfEven (Decimal p RoundHalfEven) + maclaurin x + | adjustedExponent x >= 0 = subArith (subMaclaurin x) >>= subRounded + | otherwise = sum one one one one + where sum :: FinitePrecision p + => Decimal p RoundHalfEven + -> Decimal PInfinite RoundHalfEven + -> Decimal PInfinite RoundHalfEven + -> Decimal PInfinite RoundHalfEven + -> Arith p RoundHalfEven (Decimal p RoundHalfEven) + sum s num den n = do + num' <- subArith (multiplyExact num x) + den' <- subArith (multiplyExact den n) + s' <- add s =<< divide num' den' + if s' == s then return s' + else sum s' num' den' =<< subArith (add n one) + + subMaclaurin :: FinitePrecision p => Decimal a b + -> Arith p RoundHalfEven (Decimal p RoundHalfEven) + subMaclaurin x = subArith (multiplyExact x oneHalf) >>= maclaurin >>= + \r -> multiply r r + + subRounded :: Precision p + => Decimal (PPlus1 (PPlus1 p)) a + -> Arith p r (Decimal p RoundHalfEven) + subRounded = subArith . roundDecimal + + result :: Decimal p a -> Arith p r (Decimal p a) + result r = coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r') + where r' = coerce r + +exp n@Inf { sign = s } + | s == Pos = return (coerce n) + | otherwise = return zero +exp x = coerce <$> generalRules1 x + +-- | 'fusedMultiplyAdd' takes three operands; the first two are multiplied +-- together, using 'multiply', with sufficient precision and exponent range +-- that the result is exact and unrounded. No /flags/ are set by the +-- multiplication unless one of the first two operands is a signaling NaN or +-- one is a zero and the other is an infinity. +-- +-- Unless the multiplication failed, the third operand is then added to the +-- result of that multiplication, using 'add', under the current context. +-- +-- In other words, @fusedMultiplyAdd x y z@ delivers a result which is @(x × +-- y) + z@ with only the one, final, rounding. +fusedMultiplyAdd :: (Precision p, Rounding r) + => Decimal a b -> Decimal c d -> Decimal e f + -> Arith p r (Decimal p r) +fusedMultiplyAdd x y z = + either raise (return . coerce) (exactMult x y) >>= add z + + where exactMult :: Rounding r => Decimal a b -> Decimal c d + -> Either (Exception PInfinite r) (Decimal PInfinite r) + exactMult x y = evalArith (multiply x y) newContext + + raise :: Exception a r -> Arith p r (Decimal p r) + raise e = raiseSignal (exceptionSignal e) (coerce $ exceptionResult e) + +-- | 'ln' takes one operand. If the operand is a NaN then the general rules +-- for special values apply. +-- +-- Otherwise, the operand must be a zero or positive, and the result is the +-- natural (base /e/) logarithm of the operand, with the following cases: +-- +-- * If the operand is a zero, the result is −Infinity and exact. +-- +-- * If the operand is +Infinity, the result is +Infinity and exact. +-- +-- * If the operand equals one, the result is 0 and exact. +-- +-- * Otherwise the result is inexact and will be rounded using the +-- 'RoundHalfEven' algorithm. The coefficient will have exactly /precision/ +-- digits (unless the result is subnormal). These inexact results should be +-- correctly rounded, but may be up to 1 ulp (unit in last place) in error. +ln :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p RoundHalfEven) +ln x@Num { sign = s, coefficient = c, exponent = e } + | c == 0 = return infinity { sign = Neg } + | s == Pos = if e <= 0 && c == 10^(-e) then return zero + else subArith (subLn x) >>= subRounded >>= result + + where subLn :: FinitePrecision p => Decimal a b + -> Arith p RoundHalfEven (Decimal p RoundHalfEven) + subLn x = do + let fe = fromIntegral (-(numDigits c - 1)) :: Exponent + r = fromIntegral (e - fe) :: Decimal PInfinite RoundHalfEven + lnf <- taylorLn x { exponent = fe } + add lnf =<< multiply r =<< ln10 + + subRounded :: Precision p => Decimal (PPlus1 (PPlus1 p)) a + -> Arith p r (Decimal p RoundHalfEven) + subRounded = subArith . roundDecimal + + result :: Decimal p a -> Arith p r (Decimal p a) + result r = coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r') + where r' = coerce r + +ln n@Inf { sign = Pos } = return (coerce n) +ln x = coerce <$> generalRules1 x + +taylorLn :: FinitePrecision p => Decimal a b + -> Arith p RoundHalfEven (Decimal p RoundHalfEven) +taylorLn x = do + num <- x `subtract` one + den <- x `add` one + multiply two =<< sum =<< num `divide` den + + where sum :: FinitePrecision p => Decimal p RoundHalfEven + -> Arith p RoundHalfEven (Decimal p RoundHalfEven) + sum b = multiply b b >>= \b2 -> sum' b b b2 one + + where sum' :: FinitePrecision p + => Decimal p RoundHalfEven + -> Decimal p RoundHalfEven + -> Decimal p RoundHalfEven + -> Decimal PInfinite RoundHalfEven + -> Arith p RoundHalfEven (Decimal p RoundHalfEven) + sum' s m b n = do + m' <- multiply m b + n' <- subArith (add n two) + s' <- add s =<< divide m' n' + if s' == s then return s' else sum' s' m' b n' + +ln10 :: FinitePrecision p => Arith p r (Decimal p RoundHalfEven) +ln10 = getPrecision >>= \(Just p) -> + if p <= 50 then return fastLn10 else slowLn10 + + where fastLn10 :: FinitePrecision p => Decimal p RoundHalfEven + fastLn10 = 2.3025850929940456840179914546843642076011014886288 + + slowLn10 :: FinitePrecision p => Arith p r (Decimal p RoundHalfEven) + slowLn10 = subArith (taylorLn ten) >>= subRound + + where subRound :: Precision p => Decimal (PPlus1 (PPlus1 p)) a + -> Arith p r (Decimal p RoundHalfEven) + subRound = subArith . roundDecimal + +-- | 'log10' takes one operand. If the operand is a NaN then the general rules +-- for special values apply. +-- +-- Otherwise, the operand must be a zero or positive, and the result is the +-- base 10 logarithm of the operand, with the following cases: +-- +-- * If the operand is a zero, the result is −Infinity and exact. +-- +-- * If the operand is +Infinity, the result is +Infinity and exact. +-- +-- * If the operand equals an integral power of ten (including 10^0 and +-- negative powers) and there is sufficient /precision/ to hold the integral +-- part of the result, the result is an integer (with an exponent of 0) and +-- exact. +-- +-- * Otherwise the result is inexact and will be rounded using the +-- 'RoundHalfEven' algorithm. The coefficient will have exactly /precision/ +-- digits (unless the result is subnormal). These inexact results should be +-- correctly rounded, but may be up to 1 ulp (unit in last place) in error. +log10 :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p RoundHalfEven) +log10 x@Num { sign = s, coefficient = c, exponent = e } + | c == 0 = return infinity { sign = Neg } + | s == Pos = getPrecision >>= \prec -> case powerOfTen c e of + Just p | maybe True (numDigits pc <=) prec -> return (fromInteger p) + where pc = fromInteger (Prelude.abs p) :: Coefficient + _ -> subArith (join $ divide <$> ln x <*> ln10) >>= result + + where powerOfTen :: Coefficient -> Exponent -> Maybe Integer + powerOfTen c e + | c == 10^d = Just (fromIntegral e + fromIntegral d) + | otherwise = Nothing + where d = numDigits c - 1 :: Int + + result :: Decimal p a -> Arith p r (Decimal p a) + result r = coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r') + where r' = coerce r + +log10 n@Inf { sign = Pos } = return (coerce n) +log10 x = coerce <$> generalRules1 x + +-- | 'divide' takes two operands. If either operand is a /special value/ then +-- the general rules apply. +-- +-- Otherwise, if the divisor is zero then either the Division undefined +-- condition is raised (if the dividend is zero) and the result is NaN, or the +-- Division by zero condition is raised and the result is an Infinity with a +-- sign which is the /exclusive or/ of the signs of the operands. +-- +-- Otherwise, a “long division” is effected. +-- +-- The result is then rounded to /precision/ digits, if necessary, according +-- to the /rounding/ algorithm and taking into account the remainder from the +-- division. +divide :: (FinitePrecision p, Rounding r) + => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +divide dividend@Num{ sign = xs } Num { coefficient = 0, sign = ys } + | Number.isZero dividend = divisionUndefined + | otherwise = divisionByZero infinity { sign = xorSigns xs ys } +divide Num { sign = xs, coefficient = xc, exponent = xe } + Num { sign = ys, coefficient = yc, exponent = ye } = quotient + + where quotient = result =<< answer + rn = Num { sign = rs, coefficient = rc, exponent = re } + rs = xorSigns xs ys + (rc, rem, dv, adjust) = longDivision xc yc (finitePrecision rn) + re = xe - (ye + adjust) + answer + | rem == 0 = return rn + | otherwise = roundDecimal $ case (rem * 2) `Prelude.compare` dv of + LT -> rn { coefficient = rc * 10 + 1, exponent = re - 1 } + EQ -> rn { coefficient = rc * 10 + 5, exponent = re - 1 } + GT -> rn { coefficient = rc * 10 + 9, exponent = re - 1 } + +divide Inf{} Inf{} = invalidOperation qNaN +divide Inf { sign = xs } Num { sign = ys } = + return Inf { sign = xorSigns xs ys } +divide Num { sign = xs } Inf { sign = ys } = + return zero { sign = xorSigns xs ys } +divide x y = generalRules2 x y + +type Dividend = Coefficient +type Divisor = Coefficient +type Quotient = Coefficient +type Remainder = Dividend + +longDivision :: Dividend -> Divisor -> Int + -> (Quotient, Remainder, Divisor, Exponent) +longDivision 0 dv _ = (0, 0, dv, 0) +longDivision dd dv p = step1 dd dv 0 + + where step1 :: Dividend -> Divisor -> Exponent + -> (Quotient, Remainder, Divisor, Exponent) + step1 dd dv adjust + | dd < dv = step1 (dd * 10) dv (adjust + 1) + | dd >= 10 * dv = step1 dd (dv * 10) (adjust - 1) + | otherwise = step2 dd dv adjust + + step2 :: Dividend -> Divisor -> Exponent + -> (Quotient, Remainder, Divisor, Exponent) + step2 = step3 0 + + step3 :: Quotient -> Dividend -> Divisor -> Exponent + -> (Quotient, Remainder, Divisor, Exponent) + step3 r dd dv adjust + | dv <= dd = step3 (r + 1) (dd - dv) dv adjust + | (dd == 0 && adjust >= 0) || + numDigits r == p = step4 r dd dv adjust + | otherwise = step3 (r * 10) (dd * 10) dv (adjust + 1) + + step4 :: Quotient -> Remainder -> Divisor -> Exponent + -> (Quotient, Remainder, Divisor, Exponent) + step4 = (,,,) + +reciprocal :: (FinitePrecision p, Rounding r) + => Decimal a b -> Arith p r (Decimal p r) +reciprocal = divide one + +-- | 'abs' takes one operand. If the operand is negative, the result is the +-- same as using the 'minus' operation on the operand. Otherwise, the result +-- is the same as using the 'plus' operation on the operand. +-- +-- Note that the result of this operation is affected by context and may set +-- /flags/. The 'copyAbs' operation may be used if this is not desired. +abs :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r) +abs x + | isNegative x = minus x + | otherwise = plus x + +-- | 'compare' takes two operands and compares their values numerically. If +-- either operand is a /special value/ then the general rules apply. No flags +-- are set unless an operand is a signaling NaN. +-- +-- Otherwise, the operands are compared, returning @'Right' 'LT'@ if the first +-- is less than the second, @'Right' 'EQ'@ if they are equal, or @'Right' +-- 'GT'@ if the first is greater than the second. +-- +-- A 'Left' value is returned if the result is NaN, indicating an “unordered” +-- comparison (see IEEE 754 §5.11). +compare :: Decimal a b -> Decimal c d + -> Arith p r (Either (Decimal p r) Ordering) +compare x@Num{} y@Num{} = nzp <$> subArith (subtract' xn yn) + + where subtract' :: Decimal a b -> Decimal c d + -> Arith PInfinite RoundHalfEven + (Decimal PInfinite RoundHalfEven) + subtract' = subtract + + (xn, yn) | sign x /= sign y = (either id fromOrdering $ nzp x, + either id fromOrdering $ nzp y) + | otherwise = (x, y) + + nzp :: Decimal a b -> Either (Decimal p r) Ordering + nzp Num { sign = s, coefficient = c } + | c == 0 = Right EQ + | s == Pos = Right GT + | otherwise = Right LT + nzp Inf { sign = s } = case s of + Pos -> Right GT + Neg -> Right LT + nzp n = Left (coerce n) + +compare Inf { sign = xs } Inf { sign = ys } = return $ case (xs, ys) of + (Pos, Neg) -> Right GT + (Neg, Pos) -> Right LT + _ -> Right EQ +compare Inf { sign = xs } Num { } = return $ case xs of + Pos -> Right GT + Neg -> Right LT +compare Num { } Inf { sign = ys } = return $ case ys of + Pos -> Right LT + Neg -> Right GT +compare x y = Left <$> generalRules2 x y + +-- | 'compareSignal' takes two operands and compares their values +-- numerically. This operation is identical to 'compare', except that if +-- neither operand is a signaling NaN then any quiet NaN operand is treated as +-- though it were a signaling NaN. (That is, all NaNs signal, with signaling +-- NaNs taking precedence over quiet NaNs.) +compareSignal :: Decimal a b -> Decimal c d + -> Arith p r (Either (Decimal p r) Ordering) +compareSignal x@NaN { signaling = True } y = x `compare` y +compareSignal x y@NaN { signaling = True } = x `compare` y +compareSignal x y = q2s x `compare` q2s y + + where q2s :: Decimal p r -> Decimal p r + q2s nan@NaN{} = nan { signaling = True } + q2s x = x + +-- | 'max' takes two operands, compares their values numerically, and returns +-- the maximum. If either operand is a NaN then the general rules apply, +-- unless one is a quiet NaN and the other is numeric, in which case the +-- numeric operand is returned. +max :: Decimal a b -> Decimal a b -> Arith p r (Decimal a b) +max x y = snd <$> minMax id x y + +-- | 'maxMagnitude' takes two operands and compares their values numerically +-- with their /sign/ ignored and assumed to be 0. +-- +-- If, without signs, the first operand is the larger then the original first +-- operand is returned (that is, with the original sign). If, without signs, +-- the second operand is the larger then the original second operand is +-- returned. Otherwise the result is the same as from the 'max' operation. +maxMagnitude :: Decimal a b -> Decimal a b -> Arith p r (Decimal a b) +maxMagnitude x y = snd <$> minMax withoutSign x y + +-- | 'min' takes two operands, compares their values numerically, and returns +-- the minimum. If either operand is a NaN then the general rules apply, +-- unless one is a quiet NaN and the other is numeric, in which case the +-- numeric operand is returned. +min :: Decimal a b -> Decimal a b -> Arith p r (Decimal a b) +min x y = fst <$> minMax id x y + +-- | 'minMagnitude' takes two operands and compares their values numerically +-- with their /sign/ ignored and assumed to be 0. +-- +-- If, without signs, the first operand is the smaller then the original first +-- operand is returned (that is, with the original sign). If, without signs, +-- the second operand is the smaller then the original second operand is +-- returned. Otherwise the result is the same as from the 'min' operation. +minMagnitude :: Decimal a b -> Decimal a b -> Arith p r (Decimal a b) +minMagnitude x y = fst <$> minMax withoutSign x y + +-- | Ordering function for 'min', 'minMagnitude', 'max', and 'maxMagnitude': +-- returns the original arguments as (smaller, larger) when the given function +-- is applied to them. +minMax :: (Decimal a b -> Decimal a b) -> Decimal a b -> Decimal a b + -> Arith p r (Decimal a b, Decimal a b) +minMax _ x@Num{} NaN { signaling = False } = return (x, x) +minMax _ x@Inf{} NaN { signaling = False } = return (x, x) +minMax _ NaN { signaling = False } y@Num{} = return (y, y) +minMax _ NaN { signaling = False } y@Inf{} = return (y, y) + +minMax f x y = f x `compare` f y >>= \c -> return $ case c of + Right LT -> (x, y) + Right GT -> (y, x) + Right EQ -> case (sign x, sign y) of + (Neg, Pos) -> (x, y) + (Pos, Neg) -> (y, x) + (Pos, Pos) -> case (x, y) of + (Num { exponent = xe }, Num { exponent = ye }) | xe > ye -> (y, x) + _ -> (x, y) + (Neg, Neg) -> case (x, y) of + (Num { exponent = xe }, Num { exponent = ye }) | xe < ye -> (y, x) + _ -> (x, y) + Left nan -> let nan' = coerce nan in (nan', nan') + +withoutSign :: Decimal p r -> Decimal p r +withoutSign n = n { sign = Pos } + +-- | 'power' takes two operands, and raises a number (the left-hand operand) +-- to a power (the right-hand operand). If either operand is a /special value/ +-- then the general rules apply, except as stated below. +-- +-- The following rules apply: +-- +-- * If both operands are zero, or if the left-hand operand is less than zero +-- and the right-hand operand does not have an integral value or is infinite, +-- an Invalid operation condition is raised, the result is NaN, and the +-- following rules do not apply. +-- +-- * If the left-hand operand is infinite, the result will be exact and will +-- be infinite if the right-hand side is positive, 1 if the right-hand side is +-- a zero, and 0 if the right-hand side is negative. +-- +-- * If the left-hand operand is a zero, the result will be exact and will be +-- infinite if the right-hand side is negative or 0 if the right-hand side is +-- positive. +-- +-- * If the right-hand operand is a zero, the result will be 1 and exact. +-- +-- * In cases not covered above, the result will be inexact unless the +-- right-hand side has an integral value and the result is finite and can be +-- expressed exactly within /precision/ digits. In this latter case, if the +-- result is unrounded then its exponent will be that which would result if +-- the operation were calculated by repeated multiplication (if the second +-- operand is negative then the reciprocal of the first operand is used, with +-- the absolute value of the second operand determining the multiplications). +-- +-- * Inexact finite results should be correctly rounded, but may be up to 1 +-- ulp (unit in last place) in error. +-- +-- * The /sign/ of the result will be 1 only if the right-hand side has an +-- integral value and is odd (and is not infinite) and also the /sign/ of the +-- left-hand side is 1. In all other cases, the /sign/ of the result will be +-- 0. +power :: (FinitePrecision p, Rounding r) + => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +power x@Num { coefficient = 0 } y@Num{} + | Number.isZero y = invalidOperation qNaN + | Number.isNegative y = divisionByZero infinity { sign = powerSign x y } + | otherwise = return zero { sign = powerSign x y } +power x@Num{} y@Num{} = case integralValue y of + Just i | i < 0 -> reciprocal x >>= \rx -> integralPower rx (-i) + | otherwise -> integralPower x i + Nothing | Number.isPositive x -> ln x >>= multiply y >>= fmap coerce . exp + | otherwise -> invalidOperation qNaN +power x@Num{} y@Inf{} + | Number.isPositive x = return $ case sign y of + Pos -> infinity + Neg -> zero + | otherwise = invalidOperation qNaN +power x@Inf{} y@Num{} + | Number.isZero y = return one + | Number.isPositive y = return infinity { sign = powerSign x y } + | otherwise = return zero { sign = powerSign x y } +power Inf{} Inf { sign = s } + | s == Pos = return infinity + | otherwise = return zero +power x y = generalRules2 x y + +powerSign :: Decimal a b -> Decimal c d -> Sign +powerSign x y + | Number.isNegative x && fromMaybe False (odd <$> integralValue y) = Neg + | otherwise = Pos + +integralPower :: (Precision p, Rounding r) + => Decimal a b -> Integer -> Arith p r (Decimal p r) +integralPower b e = integralPower' (return b) e one + where integralPower' :: (Precision p, Rounding r) + => Arith p r (Decimal a b) -> Integer -> Decimal p r + -> Arith p r (Decimal p r) + integralPower' _ 0 r = return r + integralPower' mb e r + | odd e = mb >>= \b -> multiply r b >>= + integralPower' (multiply b b) e' + | otherwise = integralPower' (mb >>= \b -> multiply b b) e' r + where e' = e `div` 2 + +-- | 'quantize' takes two operands. If either operand is a /special value/ +-- then the general rules apply, except that if either operand is infinite and +-- the other is finite an Invalid operation condition is raised and the result +-- is NaN, or if both are infinite then the result is the first operand. +-- +-- Otherwise (both operands are finite), 'quantize' returns the number which +-- is equal in value (except for any rounding) and sign to the first +-- (left-hand) operand and which has an /exponent/ set to be equal to the +-- exponent of the second (right-hand) operand. +-- +-- The /coefficient/ of the result is derived from that of the left-hand +-- operand. It may be rounded using the current /rounding/ setting (if the +-- /exponent/ is being increased), multiplied by a positive power of ten (if +-- the /exponent/ is being decreased), or is unchanged (if the /exponent/ is +-- already equal to that of the right-hand operand). +-- +-- Unlike other operations, if the length of the /coefficient/ after the +-- quantize operation would be greater than /precision/ then an Invalid +-- operation condition is raised. This guarantees that, unless there is an +-- error condition, the /exponent/ of the result of a quantize is always equal +-- to that of the right-hand operand. +-- +-- Also unlike other operations, quantize will never raise Underflow, even if +-- the result is subnormal and inexact. +quantize :: (Precision p, Rounding r) + => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +quantize x@Num { coefficient = xc, exponent = xe } Num { exponent = ye } + | xe > ye = result x { coefficient = xc * 10^(xe - ye), exponent = ye } + | xe < ye = rc >>= \c -> result x { coefficient = c, exponent = ye } + | otherwise = result x + + where result :: Precision p => Decimal a b -> Arith p r (Decimal p r) + result x = getPrecision >>= \p -> case numDigits (coefficient x) of + n | maybe False (n >) p -> invalidOperation qNaN + _ -> return (coerce x) + + rc :: Rounding r => Arith p r Coefficient + rc = let b = 10^(ye - xe) + (q, r) = xc `quotRem` b + in getRounder >>= \rounder -> return (rounder (sign x) r b q) + +quantize Num{} Inf{} = invalidOperation qNaN +quantize Inf{} Num{} = invalidOperation qNaN +quantize x@Inf{} Inf{} = return (coerce x) +quantize x y = generalRules2 x y + +-- | 'reduce' takes one operand. It has the same semantics as the 'plus' +-- operation, except that if the final result is finite it is reduced to its +-- simplest form, with all trailing zeros removed and its sign preserved. +reduce :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r) +reduce n = reduce' <$> plus n + where reduce' n@Num { coefficient = c, exponent = e } + | c == 0 = n { exponent = 0 } + | r == 0 = reduce' n { coefficient = q, exponent = e + 1 } + where (q, r) = c `quotRem` 10 + reduce' n = n + +-- | 'roundToIntegralExact' takes one operand. If the operand is a +-- /special value/, or the exponent of the operand is non-negative, then the +-- result is the same as the operand (unless the operand is a signaling NaN, +-- as usual). +-- +-- Otherwise (the operand has a negative exponent) the result is the same as +-- using the 'quantize' operation using the given operand as the +-- left-hand-operand, 1E+0 as the right-hand-operand, and the precision of the +-- operand as the /precision/ setting. The rounding mode is taken from the +-- context, as usual. +roundToIntegralExact :: (Precision a, Rounding r) + => Decimal a b -> Arith p r (Decimal a r) +roundToIntegralExact x@Num { exponent = e } + | e >= 0 = return (coerce x) + | otherwise = + let (Right r, context) = runArith (quantize x one) newContext + quantizeFlags = flags context + + maybeRaise :: Signal -> Decimal a r -> Arith p r (Decimal a r) + maybeRaise sig + | sig `signalMember` quantizeFlags = + fmap coerce . raiseSignal sig . coerce + | otherwise = return + + in maybeRaise Inexact r >> maybeRaise Rounded r + +roundToIntegralExact x@Inf{} = return (coerce x) +roundToIntegralExact x = coerce <$> generalRules1 x + +-- | 'roundToIntegralValue' takes one operand. It is identical to the +-- 'roundToIntegralExact' operation except that the 'Inexact' and 'Rounded' +-- flags are never set even if the operand is rounded (that is, the operation +-- is quiet unless the operand is a signaling NaN). +roundToIntegralValue :: (Precision a, Rounding r) + => Decimal a b -> Arith p r (Decimal a r) +roundToIntegralValue x@Num { exponent = e } + | e >= 0 = return (coerce x) + | otherwise = subArith (quantize x one) +roundToIntegralValue x@Inf{} = return (coerce x) +roundToIntegralValue x = coerce <$> generalRules1 x + +-- | 'squareRoot' takes one operand. If the operand is a /special value/ then +-- the general rules apply. +-- +-- Otherwise, the ideal exponent of the result is defined to be half the +-- exponent of the operand (rounded to an integer, towards −Infinity, if +-- necessary) and then: +-- +-- If the operand is less than zero an Invalid operation condition is raised. +-- +-- If the operand is greater than zero, the result is the square root of the +-- operand. If no rounding is necessary (the exact result requires /precision/ +-- digits or fewer) then the coefficient and exponent giving the correct value +-- and with the exponent closest to the ideal exponent is used. If the result +-- must be inexact, it is rounded using the 'RoundHalfEven' algorithm and the +-- coefficient will have exactly /precision/ digits (unless the result is +-- subnormal), and the exponent will be set to maintain the correct value. +-- +-- Otherwise (the operand is equal to zero), the result will be the zero with +-- the same sign as the operand and with the ideal exponent. +squareRoot :: FinitePrecision p + => Decimal a b -> Arith p r (Decimal p RoundHalfEven) +squareRoot n@Num { sign = s, coefficient = c, exponent = e } + | c == 0 = return n { exponent = idealExp } + | s == Pos = subResult >>= subRounded >>= result + + where idealExp = e `div` 2 :: Exponent + + reduced :: Decimal p r -> Decimal p r + reduced n@Num { coefficient = c, exponent = e } + | e < idealExp = case bd of + Just (b, (q, _)) -> n { coefficient = q, exponent = e + b } + Nothing -> n + | e > idealExp = n { coefficient = c * 10^d, exponent = idealExp } + where d = Prelude.abs (e - idealExp) + bd = find (\(_, (_, r)) -> r == 0) ds + ds = map (\d -> (d, c `quotRem` (10^d))) [d, d - 1 .. 1] + reduced n = n + + subResult :: FinitePrecision p + => Arith p r (Decimal (PPlus1 (PPlus1 p)) RoundHalfEven) + subResult = subArith (babylonian approx) + + subRounded :: Precision p + => Decimal a b -> Arith p r (Decimal p RoundHalfEven) + subRounded = subArith . roundDecimal + + exactness :: Decimal a b -> Arith p r + (Either (Decimal p r) Ordering) + exactness r = subArith (multiply' r r) >>= compare n + where multiply' :: Decimal a b -> Decimal c d + -> Arith PInfinite RoundHalfEven + (Decimal PInfinite RoundHalfEven) + multiply' = multiply + + result :: Decimal p a -> Arith p r (Decimal p a) + result r = exactness r >>= \e -> case e of + Right EQ -> return (reduced r) + _ -> let r' = coerce r + in coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r') + + approx :: Decimal p r + approx | even ae = n { coefficient = 2, exponent = ae `quot` 2 } + | otherwise = n { coefficient = 6, exponent = (ae - 1) `quot` 2 } + where ae = adjustedExponent n + + babylonian :: FinitePrecision p => Decimal p RoundHalfEven + -> Arith p RoundHalfEven (Decimal p RoundHalfEven) + babylonian x = do + x' <- multiply oneHalf =<< add x =<< n `divide` x + if x' == x then return x' else babylonian x' + +squareRoot n@Inf { sign = Pos } = return (coerce n) +squareRoot x = coerce <$> generalRules1 x + +-- $miscellaneous-operations +-- +-- This section describes miscellaneous operations on decimal numbers, +-- including non-numeric comparisons, sign and other manipulations, and +-- logical operations. +-- +-- Some operations return a boolean value described as 0 or 1 in the +-- /General Decimal Arithmetic Specification/, but which is returned as a +-- 'Bool' in this implementation. These values can be converted to 'Decimal' +-- via 'fromBool'. +-- +-- Similarly, the total ordering operations return an 'Ordering' value in this +-- implementation, but can be converted to 'Decimal' via 'fromOrdering'. + +data Logical = Logical { bits :: Integer, bitLength :: Int } + +toLogical :: Decimal a b -> Maybe Logical +toLogical Num { sign = Pos, coefficient = c, exponent = 0 } = + getBits c Logical { bits = zeroBits, bitLength = 0 } + + where getBits :: Coefficient -> Logical -> Maybe Logical + getBits 0 g = return g + getBits c g@Logical { bits = b, bitLength = l } = case d of + 0 -> getBits c' g { bitLength = succ l } + 1 -> getBits c' g { bits = setBit b l, bitLength = succ l } + _ -> Nothing + where (c', d) = c `quotRem` 10 + +toLogical _ = Nothing + +fromLogical :: Logical -> Decimal a b +fromLogical Logical { bits = b, bitLength = l } = + Num { sign = Pos, coefficient = fromBits 0 1 0, exponent = 0 } + + where fromBits :: Int -> Coefficient -> Coefficient -> Coefficient + fromBits i r c + | i == l = c + | testBit b i = fromBits i' r' (c + r) + | otherwise = fromBits i' r' c + where i' = succ i + r' = r * 10 + +-- $logical-operations +-- +-- The logical operations ('and', 'or', 'xor', and 'invert') take +-- /logical operands/, which are finite numbers with a /sign/ of 0, an +-- /exponent/ of 0, and a /coefficient/ whose digits must all be either 0 or +-- 1. The length of the result will be at most /precision/ digits (all of +-- which will be either 0 or 1); operands are truncated on the left or padded +-- with zeros on the left as necessary. The result of a logical operation is +-- never rounded and the only /flag/ that might be set is 'InvalidOperation' +-- (set if an operand is not a valid logical operand). + +-- | 'and' is a logical operation which takes two logical operands. The result +-- is the digit-wise /and/ of the two operands; each digit of the result is +-- the logical and of the corresponding digits of the operands, aligned at the +-- least-significant digit. A result digit is 1 if both of the corresponding +-- operand digits are 1; otherwise it is 0. +and :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +and x@Num{} y@Num{} = case (toLogical x, toLogical y) of + (Just lx, Just ly) -> getPrecision >>= \p -> + let m = Prelude.min (bitLength lx) (bitLength ly) + z = Logical { bits = bits lx .&. bits ly + , bitLength = maybe m (Prelude.min m) p } + in return (fromLogical z) + _ -> invalidOperation qNaN +and x y = generalRules2 x y + +-- | 'or' is a logical operation which takes two logical operands. The result +-- is the digit-wise /inclusive or/ of the two operands; each digit of the +-- result is the logical or of the corresponding digits of the operands, +-- aligned at the least-significant digit. A result digit is 1 if either or +-- both of the corresponding operand digits is 1; otherwise it is 0. +or :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +or x@Num{} y@Num{} = case (toLogical x, toLogical y) of + (Just lx, Just ly) -> getPrecision >>= \p -> + let m = Prelude.max (bitLength lx) (bitLength ly) + z = Logical { bits = bits lx .|. bits ly + , bitLength = maybe m (Prelude.min m) p } + in return (fromLogical z) + _ -> invalidOperation qNaN +or x y = generalRules2 x y + +-- | 'xor' is a logical operation which takes two logical operands. The result +-- is the digit-wise /exclusive or/ of the two operands; each digit of the +-- result is the logical exclusive-or of the corresponding digits of the +-- operands, aligned at the least-significant digit. A result digit is 1 if +-- one of the corresponding operand digits is 1 and the other is 0; otherwise +-- it is 0. +xor :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +xor x@Num{} y@Num{} = case (toLogical x, toLogical y) of + (Just lx, Just ly) -> getPrecision >>= \p -> + let m = Prelude.max (bitLength lx) (bitLength ly) + z = Logical { bits = bits lx `Bits.xor` bits ly + , bitLength = maybe m (Prelude.min m) p } + in return (fromLogical z) + _ -> invalidOperation qNaN +xor x y = generalRules2 x y + +-- | 'invert' is a logical operation which takes one logical operand. The +-- result is the digit-wise /inversion/ of the operand; each digit of the +-- result is the inverse of the corresponding digit of the operand. A result +-- digit is 1 if the corresponding operand digit is 0; otherwise it is 0. +invert :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p r) +invert x@Num{} = case toLogical x of + Just lx -> getPrecision >>= \(Just p) -> + let z = Logical { bits = complement (bits lx), bitLength = p } + in return (fromLogical z) + _ -> invalidOperation qNaN +invert x = generalRules1 x + +-- | 'canonical' takes one operand. The result has the same value as the +-- operand but always uses a /canonical/ encoding. The definition of +-- /canonical/ is implementation-defined; if more than one internal encoding +-- for a given NaN, Infinity, or finite number is possible then one +-- “preferred” encoding is deemed canonical. This operation then returns the +-- value using that preferred encoding. +-- +-- If all possible operands have just one internal encoding each, then +-- 'canonical' always returns the operand unchanged (that is, it has the same +-- effect as 'copy'). This operation is unaffected by context and is quiet — +-- no /flags/ are changed in the context. +canonical :: Decimal a b -> Arith p r (Decimal a b) +canonical = copy + +-- | 'class'' takes one operand. The result is an indication of the /class/ of +-- the operand, where the class is one of ten possibilities, corresponding to +-- one of the strings @"sNaN"@ (signaling NaN), @\"NaN"@ (quiet NaN), +-- @"-Infinity"@ (negative infinity), @"-Normal"@ (negative normal finite +-- number), @"-Subnormal"@ (negative subnormal finite number), @"-Zero"@ +-- (negative zero), @"+Zero"@ (non-negative zero), @"+Subnormal"@ (positive +-- subnormal finite number), @"+Normal"@ (positive normal finite number), or +-- @"+Infinity"@ (positive infinity). This operation is quiet; no /flags/ are +-- changed in the context. +-- +-- Note that unlike the special values in the model, the sign of any NaN is +-- ignored in the classification, as required by IEEE 754. +class' :: Precision a => Decimal a b -> Arith p r Class +class' n = return $ case n of + Num {} | Number.isZero n -> NumberClass (sign n) ZeroClass + | Number.isSubnormal n -> NumberClass (sign n) SubnormalClass + | otherwise -> NumberClass (sign n) NormalClass + Inf {} -> NumberClass (sign n) InfinityClass + NaN { signaling = s } -> NaNClass (toEnum . fromEnum $ s) + +data Class = NumberClass Sign NumberClass -- ^ Number (finite or infinite) + | NaNClass NaNClass -- ^ Not a number (quiet or signaling) + deriving Eq + +data NumberClass = ZeroClass -- ^ Zero + | SubnormalClass -- ^ Subnormal finite number + | NormalClass -- ^ Normal finite number + | InfinityClass -- ^ Infinity + deriving Eq + +data NaNClass = QuietClass -- ^ Quiet NaN + | SignalingClass -- ^ Signaling NaN + deriving (Eq, Enum) + +instance Show Class where + show c = case c of + NumberClass s nc -> signChar s : showNumberClass nc + NaNClass QuietClass -> nan + NaNClass SignalingClass -> 's' : nan + + where signChar :: Sign -> Char + signChar Pos = '+' + signChar Neg = '-' + + showNumberClass :: NumberClass -> String + showNumberClass nc = case nc of + ZeroClass -> "Zero" + SubnormalClass -> "Subnormal" + NormalClass -> "Normal" + InfinityClass -> "Infinity" + + nan :: String + nan = "NaN" + +-- | 'compareTotal' takes two operands and compares them using their abstract +-- representation rather than their numerical value. A /total ordering/ is +-- defined for all possible abstract representations, as described below. If +-- the first operand is lower in the total order than the second operand then +-- the result is 'LT', if the operands have the same abstract representation +-- then the result is 'EQ', and if the first operand is higher in the total +-- order than the second operand then the result is 'GT'. The total ordering +-- is defined as follows. +-- +-- 1. The following items describe the ordering for representations whose +-- /sign/ is 0. If the /sign/ is 1, the order is reversed. A representation +-- with a /sign/ of 1 is always lower in the ordering than one with a /sign/ +-- of 0. +-- +-- 2. Numbers (representations which are not NaNs) are ordered such that a +-- larger numerical value is higher in the ordering. If two representations +-- have the same numerical value then the exponent is taken into account; +-- larger (more positive) exponents are higher in the ordering. +-- +-- 3. All quiet NaNs are higher in the total ordering than all signaling NaNs. +-- +-- 4. Quiet NaNs and signaling NaNs are ordered according to their /payload/; +-- a larger payload is higher in the ordering. +-- +-- For example, the following values are ordered from lowest to highest: @-NaN +-- -sNaN -Infinity -127 -1 -1.00 -0 -0.000 0 1.2300 1.23 1E+9 Infinity sNaN +-- NaN NaN456@. +compareTotal :: Decimal a b -> Decimal c d -> Arith p r Ordering +compareTotal x y = return $ case (sign x, sign y) of + (Pos, Pos) -> compareAbs x y + (Neg, Neg) -> compareAbs y x + (Neg, Pos) -> LT + (Pos, Neg) -> GT + + where compareAbs :: Decimal a b -> Decimal c d -> Ordering + compareAbs Num { coefficient = xc, exponent = xe } + Num { coefficient = yc, exponent = ye } = + let (xac, yac) | xe == ye = (xc, yc) + | xe > ye = (xc * 10^n, yc) + | otherwise = (xc, yc * 10^n) + n = Prelude.abs (xe - ye) + in Prelude.compare xac yac `mappend` Prelude.compare xe ye + compareAbs Num{} Inf{} = LT + compareAbs Inf{} Num{} = GT + compareAbs Inf{} Inf{} = EQ + compareAbs NaN { signaling = xs, payload = xp } + NaN { signaling = ys, payload = yp } = + Prelude.compare ys xs `mappend` Prelude.compare xp yp + compareAbs NaN{} _ = GT + compareAbs _ NaN{} = LT + +-- | 'compareTotalMagnitude' takes two operands and compares them using their +-- abstract representation rather than their numerical value and with their +-- /sign/ ignored and assumed to be 0. The result is identical to that +-- obtained by using 'compareTotal' on two operands which are the 'copyAbs' +-- copies of the operands to 'compareTotalMagnitude'. +compareTotalMagnitude :: Decimal a b -> Decimal c d -> Arith p r Ordering +compareTotalMagnitude x y = compareTotal x { sign = Pos } y { sign = Pos } + +-- | 'copy' takes one operand. The result is a copy of the operand. This +-- operation is unaffected by context and is quiet — no /flags/ are changed in +-- the context. +copy :: Decimal a b -> Arith p r (Decimal a b) +copy = return + +-- | 'copyAbs' takes one operand. The result is a copy of the operand with the +-- /sign/ set to 0. Unlike the 'abs' operation, this operation is unaffected +-- by context and is quiet — no /flags/ are changed in the context. +copyAbs :: Decimal a b -> Arith p r (Decimal a b) +copyAbs n = return n { sign = Pos } + +-- | 'copyNegate' takes one operand. The result is a copy of the operand with +-- the /sign/ inverted (a /sign/ of 0 becomes 1 and vice versa). Unlike the +-- 'minus' operation, this operation is unaffected by context and is quiet — +-- no /flags/ are changed in the context. +copyNegate :: Decimal a b -> Arith p r (Decimal a b) +copyNegate n = return n { sign = negateSign (sign n) } + +-- | 'copySign' takes two operands. The result is a copy of the first operand +-- with the /sign/ set to be the same as the /sign/ of the second +-- operand. This operation is unaffected by context and is quiet — no /flags/ +-- are changed in the context. +copySign :: Decimal a b -> Decimal c d -> Arith p r (Decimal a b) +copySign n m = return n { sign = sign m } + +-- | 'isCanonical' takes one operand. The result is 'True' if the operand is +-- /canonical/; otherwise it is 'False'. The definition of /canonical/ is +-- implementation-defined; if more than one internal encoding for a given NaN, +-- Infinity, or finite number is possible then one “preferred” encoding is +-- deemed canonical. This operation then tests whether the internal encoding +-- is that preferred encoding. +-- +-- If all possible operands have just one internal encoding each, then +-- 'isCanonical' always returns 'True'. This operation is unaffected by +-- context and is quiet — no /flags/ are changed in the context. +isCanonical :: Decimal a b -> Arith p r Bool +isCanonical _ = return True + +-- | 'isFinite' takes one operand. The result is 'True' if the operand is +-- neither infinite nor a NaN (that is, it is a normal number, a subnormal +-- number, or a zero); otherwise it is 'False'. This operation is unaffected +-- by context and is quiet — no /flags/ are changed in the context. +isFinite :: Decimal a b -> Arith p r Bool +isFinite = return . Number.isFinite + +-- | 'isInfinite' takes one operand. The result is 'True' if the operand is an +-- Infinity; otherwise it is 'False'. This operation is unaffected by context +-- and is quiet — no /flags/ are changed in the context. +isInfinite :: Decimal a b -> Arith p r Bool +isInfinite n = return $ case n of + Inf{} -> True + _ -> False + +-- | 'isNaN' takes one operand. The result is 'True' if the operand is a NaN +-- (quiet or signaling); otherwise it is 'False'. This operation is unaffected +-- by context and is quiet — no /flags/ are changed in the context. +isNaN :: Decimal a b -> Arith p r Bool +isNaN n = return $ case n of + NaN{} -> True + _ -> False + +-- | 'isNormal' takes one operand. The result is 'True' if the operand is a +-- positive or negative /normal number/; otherwise it is 'False'. This +-- operation is quiet; no /flags/ are changed in the context. +isNormal :: Precision a => Decimal a b -> Arith p r Bool +isNormal = return . Number.isNormal + +-- | 'isQNaN' takes one operand. The result is 'True' if the operand is a +-- quiet NaN; otherwise it is 'False'. This operation is unaffected by context +-- and is quiet — no /flags/ are changed in the context. +isQNaN :: Decimal a b -> Arith p r Bool +isQNaN n = return $ case n of + NaN { signaling = False } -> True + _ -> False + +-- | 'isSigned' takes one operand. The result is 'True' if the /sign/ of the +-- operand is 1; otherwise it is 'False'. This operation is unaffected by +-- context and is quiet — no /flags/ are changed in the context. +isSigned :: Decimal a b -> Arith p r Bool +isSigned = return . Number.isNegative + +-- | 'isSNaN' takes one operand. The result is 'True' if the operand is a +-- signaling NaN; otherwise it is 'False'. This operation is unaffected by +-- context and is quiet — no /flags/ are changed in the context. +isSNaN :: Decimal a b -> Arith p r Bool +isSNaN n = return $ case n of + NaN { signaling = True } -> True + _ -> False + +-- | 'isSubnormal' takes one operand. The result is 'True' if the operand is a +-- positive or negative /subnormal number/; otherwise it is 'False'. This +-- operation is quiet; no /flags/ are changed in the context. +isSubnormal :: Precision a => Decimal a b -> Arith p r Bool +isSubnormal = return . Number.isSubnormal + +-- | 'isZero' takes one operand. The result is 'True' if the operand is a +-- zero; otherwise it is 'False'. This operation is unaffected by context and +-- is quiet — no /flags/ are changed in the context. +isZero :: Decimal a b -> Arith p r Bool +isZero = return . Number.isZero + +-- | 'logb' takes one operand. If the operand is a NaN then the general +-- arithmetic rules apply. If the operand is infinite then +Infinity is +-- returned. If the operand is a zero, then −Infinity is returned and the +-- Division by zero exceptional condition is raised. +-- +-- Otherwise, the result is the integer which is the exponent of the magnitude +-- of the most significant digit of the operand (as though the operand were +-- truncated to a single digit while maintaining the value of that digit and +-- without limiting the resulting exponent). All results are exact unless an +-- integer result does not fit in the available /precision/. +logb :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r) +logb Num { coefficient = c, exponent = e } + | c == 0 = raiseSignal DivisionByZero Inf { sign = Neg } + | otherwise = roundDecimal (fromInteger r :: Decimal PInfinite RoundHalfEven) + where r = fromIntegral (numDigits c) - 1 + fromIntegral e :: Integer +logb Inf{} = return Inf { sign = Pos } +logb x = generalRules1 x + +-- | 'scaleb' takes two operands. If either operand is a NaN then the general +-- arithmetic rules apply. Otherwise, the second operand must be a finite +-- integer with an exponent of zero and in the range ±2 × (E/max/ + +-- /precision/) inclusive, where E/max/ is the largest value that can be +-- returned by the 'logb' operation at the same /precision/ setting. (If is is +-- not, the Invalid Operation condition is raised and the result is NaN.) +-- +-- If the first operand is infinite then that Infinity is returned, otherwise +-- the result is the first operand modified by adding the value of the second +-- operand to its /exponent/. The result may Overflow or Underflow. +scaleb :: Decimal a b -> Decimal c d -> Arith p r (Decimal a b) +scaleb x@Num { exponent = e } s + | validScale s = let Just i = integralValue s + in return x { exponent = e + fromInteger i } + -- XXX check for Overflow and Underflow +scaleb x@Inf{} s | validScale s = return x +scaleb x y = coerce <$> generalRules2 x y + +validScale :: Decimal a b -> Bool +validScale Num { exponent = 0 } = True -- XXX +validScale _ = False + +-- | 'radix' takes no operands. The result is the radix (base) in which +-- arithmetic is effected; for this specification the result will have the +-- value 10. +radix :: Precision p => Arith p r (Decimal p r) +radix = return radix' + where radix' = case precision radix' of + Just 1 -> one { exponent = 1 } + _ -> one { coefficient = 10 } + +-- | 'sameQuantum' takes two operands, and returns 'True' if the two operands +-- have the same /exponent/ or 'False' otherwise. The result is never affected +-- by either the sign or the coefficient of either operand. +-- +-- If either operand is a /special value/, 'True' is returned only if both +-- operands are NaNs or both are infinities. +-- +-- 'sameQuantum' does not change any /flags/ in the context. +sameQuantum :: Decimal a b -> Decimal c d -> Arith p r Bool +sameQuantum Num { exponent = xe } Num { exponent = ye } = return (xe == ye) +sameQuantum Inf { } Inf { } = return True +sameQuantum NaN { } NaN { } = return True +sameQuantum _ _ = return False + +-- | 'shift' takes two operands. The second operand must be an integer (with +-- an /exponent/ of 0) in the range /−precision/ through /precision/. If the +-- first operand is a NaN then the general arithmetic rules apply, and if it +-- is infinite then the result is the Infinity unchanged. +-- +-- Otherwise (the first operand is finite) the result has the same /sign/ and +-- /exponent/ as the first operand, and a /coefficient/ which is a shifted +-- copy of the digits in the coefficient of the first operand. The number of +-- places to shift is taken from the absolute value of the second operand, +-- with the shift being to the left if the second operand is positive or to +-- the right otherwise. Digits shifted into the coefficient are zeros. +-- +-- The only /flag/ that might be set is 'InvalidOperation' (set if the first +-- operand is an sNaN or the second is not valid). +-- +-- The 'rotate' operation can be used to rotate rather than shift a +-- coefficient. +shift :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +shift n@Num { coefficient = c } s@Num { sign = d, coefficient = sc } + | validShift z s = return z + where z = case precision z of + Just p -> y { coefficient = coefficient y `rem` 10 ^ p } + Nothing -> y + y = case d of + Pos -> n { coefficient = c * 10 ^ sc } + Neg -> n { coefficient = c `quot` 10 ^ sc } + +shift n@Inf { } s | validShift z s = return z + where z = coerce n +shift n@NaN { signaling = False } s | validShift z s = return z + where z = coerce n +shift n@NaN { signaling = True } _ = invalidOperation n +shift _ s = invalidOperation s + +validShift :: Precision p => p -> Decimal a b -> Bool +validShift px Num { coefficient = c, exponent = 0 } = + let p = fromIntegral <$> precision px in maybe True (c <=) p +validShift _ _ = False + +-- | 'rotate' takes two operands. The second operand must be an integer (with +-- an /exponent/ of 0) in the range /−precision/ through /precision/. If the +-- first operand is a NaN then the general arithmetic rules apply, and if it +-- is infinite then the result is the Infinity unchanged. +-- +-- Otherwise (the first operand is finite) the result has the same /sign/ and +-- /exponent/ as the first operand, and a /coefficient/ which is a rotated +-- copy of the digits in the coefficient of the first operand. The number of +-- places of rotation is taken from the absolute value of the second operand, +-- with the rotation being to the left if the second operand is positive or to +-- the right otherwise. +-- +-- If the coefficient of the first operand has fewer than /precision/ digits, +-- it is treated as though it were padded on the left with zeros to length +-- /precision/ before the rotation. Similarly, if the coefficient of the first +-- operand has more than /precision/ digits, it is truncated on the left +-- before use. +-- +-- The only /flag/ that might be set is 'InvalidOperation' (set if the first +-- operand is an sNaN or the second is not valid). +-- +-- The 'shift' operation can be used to shift rather than rotate a +-- coefficient. +rotate :: FinitePrecision p + => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +rotate n@Num { coefficient = c } s@Num { sign = d, coefficient = sc } + | validShift z s = return z + where z = n { coefficient = rc * b + (lc `rem` b) } + (lc, rc) = c `quotRem` b' + (b , b') = case d of + Pos -> (10^sc , 10^sc') + Neg -> (10^sc', 10^sc ) + sc' = finitePrecision z - fromIntegral sc + +rotate n@Inf { } s | validShift z s = return z + where z = coerce n +rotate n@NaN { signaling = False } s | validShift z s = return z + where z = coerce n +rotate n@NaN { signaling = True } _ = invalidOperation n +rotate _ s = invalidOperation s diff --git a/src/Numeric/Decimal/Operation.hs-boot b/src/Numeric/Decimal/Operation.hs-boot --- a/src/Numeric/Decimal/Operation.hs-boot +++ b/src/Numeric/Decimal/Operation.hs-boot @@ -11,6 +11,7 @@ , minus , abs , compare + , compareTotal , min , max , power @@ -48,12 +49,11 @@ => Decimal a b -> Arith p r (Decimal p r) abs :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r) -compare :: (Precision p, Rounding r) - => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) -min :: (Precision p, Rounding r) - => Decimal a b -> Decimal a b -> Arith p r (Decimal a b) -max :: (Precision p, Rounding r) - => Decimal a b -> Decimal a b -> Arith p r (Decimal a b) +compare :: Decimal a b -> Decimal c d + -> Arith p r (Either (Decimal p r) Ordering) +compareTotal :: Decimal a b -> Decimal c d -> Arith p r Ordering +min :: Decimal a b -> Decimal a b -> Arith p r (Decimal a b) +max :: Decimal a b -> Decimal a b -> Arith p r (Decimal a b) power :: (FinitePrecision p, Rounding r) => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) squareRoot :: FinitePrecision p diff --git a/src/Numeric/Decimal/Precision.hs b/src/Numeric/Decimal/Precision.hs --- a/src/Numeric/Decimal/Precision.hs +++ b/src/Numeric/Decimal/Precision.hs @@ -16,12 +16,27 @@ , PInfinite ) where +import {-# SOURCE #-} Numeric.Decimal.Number + -- | Precision indicates the maximum number of significant decimal digits a -- number may have. class Precision p where - -- | Return the precision of the argument, or 'Nothing' if the precision is infinite. + -- | Return the precision of the argument, or 'Nothing' if the precision is + -- infinite. precision :: p -> Maybe Int + -- | Return the maximum exponent for a number in scientific notation with + -- the given precision, or 'Nothing' if the exponent has no limit. + eMax :: p -> Maybe Exponent + eMax n = subtract 1 . (10 ^) . numDigits <$> base + where mlength = precision n :: Maybe Int + base = (10 *) . fromIntegral <$> mlength :: Maybe Coefficient + + -- | Return the minimum exponent for a number in scientific notation with + -- the given precision, or 'Nothing' if the exponent has no limit. + eMin :: p -> Maybe Exponent + eMin = fmap (1 -) . eMax + -- | A subclass of precisions that are finite class Precision p => FinitePrecision p @@ -29,6 +44,7 @@ data PInfinite instance Precision PInfinite where precision _ = Nothing + eMax _ = Nothing -- | A precision of 1 significant digit data P1 diff --git a/src/Numeric/Decimal/Rounding.hs b/src/Numeric/Decimal/Rounding.hs --- a/src/Numeric/Decimal/Rounding.hs +++ b/src/Numeric/Decimal/Rounding.hs @@ -1,7 +1,7 @@ module Numeric.Decimal.Rounding ( RoundingAlgorithm(..) - , Rounding(..) + , Rounding(rounding) , RoundDown , RoundHalfUp @@ -14,6 +14,8 @@ , Round05Up , getRounder + , Rounder + , roundDecimal ) where @@ -21,9 +23,10 @@ import Data.Coerce (coerce) +import {-# SOURCE #-} Numeric.Decimal.Arithmetic +import {-# SOURCE #-} Numeric.Decimal.Exception import {-# SOURCE #-} Numeric.Decimal.Number import Numeric.Decimal.Precision -import {-# SOURCE #-} Numeric.Decimal.Arithmetic -- | A value representation of a rounding algorithm (cf. 'Rounding'). data RoundingAlgorithm = RoundDown @@ -53,7 +56,7 @@ getRounder' = return roundCoefficient -- | Round a 'Decimal' to the precision of the arithmetic context using the --- rounding algorithm of the arithmetic context. +-- rounding mode of the arithmetic context. roundDecimal :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r) roundDecimal n@Num { sign = s, coefficient = c, exponent = e } = do @@ -68,20 +71,23 @@ n' = case excessDigits c' =<< p of Nothing -> n { coefficient = c' , exponent = e' } _ -> n { coefficient = c' `quot` 10, exponent = succ e' } - rounded :: Decimal p r -> Arith p r (Decimal p r) - rounded - | r /= 0 = raiseSignal Inexact - | otherwise = return - raiseSignal Rounded =<< rounded n' -- XXX check for overflow + rounded =<< (if r /= 0 then inexact else return) n' + -- XXX check for overflow Nothing -> return (coerce n) - where excessDigits :: Coefficient -> Int -> Maybe Int - excessDigits c p | d > p = Just (d - p) - | otherwise = Nothing - where d = numDigits c :: Int +roundDecimal n@NaN { payload = p } = do + prec <- getPrecision + case excessDigits p =<< (pred <$> prec) of + Just _ -> return n { payload = 0 } + Nothing -> return (coerce n) roundDecimal n = return (coerce n) + +excessDigits :: Coefficient -> Int -> Maybe Int +excessDigits c p | d > p = Just (d - p) + | otherwise = Nothing + where d = numDigits c :: Int -- Required algorithms... diff --git a/stack.yaml b/stack.yaml --- a/stack.yaml +++ b/stack.yaml @@ -1,13 +1,13 @@ # This file was automatically generated by 'stack init' -# +# # Some commonly used options have been documented as comments in this file. # For advanced use and comprehensive documentation of the format, please see: -# http://docs.haskellstack.org/en/stable/yaml_configuration/ +# https://docs.haskellstack.org/en/stable/yaml_configuration/ # Resolver to choose a 'specific' stackage snapshot or a compiler version. # A snapshot resolver dictates the compiler version and the set of packages # to be used for project dependencies. For example: -# +# # resolver: lts-3.5 # resolver: nightly-2015-09-21 # resolver: ghc-7.10.2 @@ -15,11 +15,11 @@ # resolver: # name: custom-snapshot # location: "./custom-snapshot.yaml" -resolver: lts-6.9 +resolver: lts-9.1 # User packages to be built. # Various formats can be used as shown in the example below. -# +# # packages: # - some-directory # - https://example.com/foo/bar/baz-0.0.2.tar.gz @@ -31,12 +31,12 @@ # subdirs: # - auto-update # - wai -# +# # A package marked 'extra-dep: true' will only be built if demanded by a # non-dependency (i.e. a user package), and its test suites and benchmarks # will not be run. This is useful for tweaking upstream packages. packages: -- '.' +- . # Dependency packages to be pulled from upstream that are not in the resolver # (e.g., acme-missiles-0.3) extra-deps: [] @@ -49,18 +49,18 @@ # Control whether we use the GHC we find on the path # system-ghc: true -# +# # Require a specific version of stack, using version ranges # require-stack-version: -any # Default -# require-stack-version: ">=1.1" -# +# require-stack-version: ">=1.5" +# # Override the architecture used by stack, especially useful on Windows # arch: i386 # arch: x86_64 -# +# # Extra directories used by stack for building # extra-include-dirs: [/path/to/dir] # extra-lib-dirs: [/path/to/dir] -# +# # Allow a newer minor version of GHC than the snapshot specifies # compiler-check: newer-minor diff --git a/test/Arbitrary.hs b/test/Arbitrary.hs new file mode 100644 --- /dev/null +++ b/test/Arbitrary.hs @@ -0,0 +1,31 @@ + +module Arbitrary + ( Arbitrary + ) where + +import Numeric.Decimal +import Test.QuickCheck + +infinity :: (Precision p, Rounding r) => Decimal p r +infinity = read "Infinity" + +instance (Precision p, Rounding r) => Arbitrary (Decimal p r) where + arbitrary = frequency [(85, genNum), (10, genInf)] + +genNum :: (Precision p, Rounding r) => Gen (Decimal p r) +genNum = do + c <- choose (-(10^10), 10^10) :: Gen Integer + e <- choose (-99, 99) :: Gen Integer + return $ read (show c ++ 'E' : show e) + +genInf :: (Precision p, Rounding r) => Gen (Decimal p r) +genInf = do + s <- elements [-1, 1] + return (s * infinity) + +genNaN :: (Precision p, Rounding r) => Gen (Decimal p r) +genNaN = oneof [nan "", nan "s"] + where nan kind = do + s <- elements ["", "-"] + p <- choose (0, 10000) :: Gen Integer + return $ read (s ++ kind ++ "NaN" ++ show p) diff --git a/test/Numeric/Decimal/EncodingSpec.hs b/test/Numeric/Decimal/EncodingSpec.hs new file mode 100644 --- /dev/null +++ b/test/Numeric/Decimal/EncodingSpec.hs @@ -0,0 +1,25 @@ + +{-# LANGUAGE OverloadedStrings #-} + +module Numeric.Decimal.EncodingSpec (spec) where + +import Test.Hspec +import Data.Binary + +import Numeric.Decimal.Encoding + +spec :: Spec +spec = do + it "encodes (-7.50) correctly" $ + encode (read "-7.50" :: Decimal64) `shouldBe` + "\xA2\x30\x00\x00\x00\x00\x03\xD0" + it "decodes (-7.50) correctly" $ + (decode "\xA2\x30\x00\x00\x00\x00\x03\xD0" :: Decimal64) `shouldBe` (-7.50) + + it "decodes Infinity correctly" $ + (decode "\x78\xFF\xFF\xFF\xFF\xFF\xFF\xFF" :: Decimal64) `shouldSatisfy` + \x -> isInfinite x && signum x == 1 + +-- prop> decode (encode x) == (x :: Decimal32) +-- prop> decode (encode x) == (x :: Decimal64) +-- prop> decode (encode x) == (x :: Decimal128) diff --git a/test/Numeric/Decimal/NumberSpec.hs b/test/Numeric/Decimal/NumberSpec.hs new file mode 100644 --- /dev/null +++ b/test/Numeric/Decimal/NumberSpec.hs @@ -0,0 +1,165 @@ + +module Numeric.Decimal.NumberSpec (spec) where + +import Arbitrary () + +import Test.Hspec +import Test.QuickCheck +import Numeric.Decimal + +spec :: Spec +spec = do + describe "Read" $ do + it "reads with Just correctly (positive)" $ + (read "Just 123" :: Maybe GeneralDecimal) `shouldBe` Just 123 + it "reads with Just correctly (negative)" $ + (read "Just (-12.0)" :: Maybe GeneralDecimal) `shouldBe` Just (-12) + + describe "Ord" $ do + it "satisfies (>) invariant" $ + property $ \x y -> + x > y ==> max x y == x && max y x == (x :: BasicDecimal) + it "satisfies (<) invariant" $ + property $ \x y -> + x < y ==> min x y == x && min y x == (x :: BasicDecimal) + + it "satisfies `max` invariant with respect to 1st arg" $ + property $ \x y -> + max x y == x ==> x >= (y :: BasicDecimal) + it "satisfies `max` invariant with respect to 2nd arg" $ + property $ \x y -> + max x y == y ==> y >= (x :: BasicDecimal) + it "satisfies `min` invariant with respect to 1st arg" $ + property $ \x y -> + min x y == x ==> x <= (y :: BasicDecimal) + it "satisfies `min` invariant with respect to 2nd arg" $ + property $ \x y -> + min x y == y ==> y <= (x :: BasicDecimal) + + describe "Enum" $ do + it "enumerates `enumFromTo` precisely" $ + ([1.7 .. 5.7] :: [BasicDecimal]) `shouldBe` [1.7, 2.7, 3.7, 4.7, 5.7] + + it "enumerates `enumFromThenTo` precisely (ascending)" $ + ([0, 0.1 .. 2] :: [BasicDecimal]) `shouldBe` + [ 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, + 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0] + it "enumerates `enumFromThenTo` precisely (descending)" $ + ([2, 1.9 .. 0] :: [BasicDecimal]) `shouldBe` + [ 2, 1.9, 1.8, 1.7, 1.6, 1.5, 1.4, 1.3, 1.2, 1.1, + 1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.0] + + describe "Num" $ do + it "satisfies relation between (+) and (*)" $ + property $ \x -> + x + x == x * (2 :: GeneralDecimal) + it "satisfies relation between (-) and 0" $ + property $ \x -> + isFinite x ==> x - x == (0 :: BasicDecimal) + it "satisfies relation between (+) and `negate`" $ + property $ \x -> + isFinite x ==> x + negate x == (0 :: BasicDecimal) + it "satisfies relation between `abs` and 0" $ + property $ \x -> + abs x >= (0 :: GeneralDecimal) + it "satisfies relation between `abs` and `signum`" $ + property $ \x -> + abs x * signum x == (x :: GeneralDecimal) + + describe "Fractional" $ do + context "RoundHalfUp" $ do + it "properly rounds (down)" $ + (4.14 :: Decimal P2 RoundHalfUp) `shouldBe` 4.1 + it "properly rounds (up)" $ + (4.15 :: Decimal P2 RoundHalfUp) `shouldBe` 4.2 + context "RoundHalfDown" $ + it "properly rounds" $ + (4.15 :: Decimal P2 RoundHalfDown) `shouldBe` 4.1 + context "RoundHalfEven" $ do + it "properly rounds 1 (up)" $ + (4.15 :: Decimal P2 RoundHalfEven) `shouldBe` 4.2 + it "properly rounds 1 (down)" $ + (4.25 :: Decimal P2 RoundHalfEven) `shouldBe` 4.2 + it "properly rounds 2 (up)" $ + (4.35 :: Decimal P2 RoundHalfEven) `shouldBe` 4.4 + it "properly rounds 2 (down)" $ + (4.45 :: Decimal P2 RoundHalfEven) `shouldBe` 4.4 + + describe "RealFrac" $ do + it "satisfies `properFraction` invariant 1" $ + property $ \x -> + let (n,f) = properFraction (x :: BasicDecimal) :: (Integer, BasicDecimal) + in x == fromIntegral n + f + it "satisfies `properFraction` invariant 2" $ + property $ \x -> + let (n,_) = properFraction (x :: BasicDecimal) :: (Integer, BasicDecimal) + in (x < 0 && n <= 0) || (x >= 0 && n >= 0) + it "satisfies `properFraction` invariant 3" $ + property $ \x -> + let (_,f) = properFraction (x :: BasicDecimal) :: (Integer, BasicDecimal) + in (x < 0 && f <= 0) || (x >= 0 && f >= 0) + it "satisfies `properFraction` invariant 4" $ + property $ \x -> + let (_,f) = properFraction (x :: BasicDecimal) :: (Integer, BasicDecimal) + in isFinite f ==> abs f < 1 + + describe "Floating" $ do + it "produces same `pi` as Double" $ + realToFrac (pi :: ExtendedDecimal P16) `shouldBe` (pi :: Double) + it "satisfies relation between (**) and (^)" $ + property $ \x y -> + y >= 0 ==> (x :: BasicDecimal) ** fromInteger y == x ^ y + it "computes `sqrt` correctly" $ + property $ \x -> + isFinite x && x >= 0 ==> + sqrt (x * x) `shouldBe` (x :: ExtendedDecimal P32) + -- coefficient (sqrt (x * x) - (x :: ExtendedDecimal P16)) <= 1 + + describe "RealFloat" $ do + it "satisfies `decodeFloat` invariant 1" $ + property $ \x -> + isFinite x ==> let b = floatRadix (x :: BasicDecimal) + (m, n) = decodeFloat x + in x == fromInteger m * fromInteger b ^^ n + it "satisfies `decodeFloat` invariant 2 (zero)" $ + decodeFloat (0 :: BasicDecimal) `shouldBe` (0,0) + it "satisfies `decodeFloat` invariant 2 (negative zero)" $ + decodeFloat (read "-0" :: BasicDecimal) `shouldBe` (0,0) + it "satisfies `decodeFloat` invariant 2 (nonzero)" $ + property $ \x -> + isFinite x && x /= 0 ==> let b = floatRadix (x :: BasicDecimal) + (m, _) = decodeFloat x + d = floatDigits x + am = abs m + in b^(d-1) <= am && am < b^d + it "satisfies relation between `encodeFloat` and `decodeFloat`" $ + property $ \x -> + not (isNegativeZero x) ==> + uncurry encodeFloat (decodeFloat x) == (x :: BasicDecimal) + it "satisfies `exponent` invariant 1 (zero)" $ + exponent (0 :: BasicDecimal) `shouldBe` 0 + it "satisfies `exponent` invariant 1 (nonzero)" $ + property $ \x -> + isFinite x && x /= 0 ==> + exponent (x :: BasicDecimal) == snd (decodeFloat x) + floatDigits x + it "satisfies `exponent` invariant 2" $ + property $ \x -> + isFinite x ==> let b = floatRadix (x :: BasicDecimal) + in x == significand x * fromInteger b ^^ exponent x + it "satisfies `significand` invariant" $ + property $ \x -> + isFinite x ==> let s = significand (x :: BasicDecimal) + b = floatRadix x + in s == 0 || + (s > -1 && s < 1 && abs s >= 1 / fromInteger b) + + it "detects negative zero" $ + isNegativeZero (read "-0" :: BasicDecimal) `shouldBe` True + it "does not detect normal zero as negative" $ + isNegativeZero (read "+0" :: BasicDecimal) `shouldBe` False + it "does not detect nonzero numbers as negative zero" $ + property $ \x -> + x /= 0 ==> isNegativeZero (x :: BasicDecimal) == False + +isFinite :: (FinitePrecision p, Rounding r) => Decimal p r -> Bool +isFinite x = not (isNaN x || isInfinite x) diff --git a/test/Numeric/Decimal/OperationSpec.hs b/test/Numeric/Decimal/OperationSpec.hs new file mode 100644 --- /dev/null +++ b/test/Numeric/Decimal/OperationSpec.hs @@ -0,0 +1,618 @@ + +module Numeric.Decimal.OperationSpec (spec) where + +import Test.Hspec + +import Numeric.Decimal +import Numeric.Decimal.Arithmetic + +import qualified Numeric.Decimal.Operation as Op + +spec :: Spec +spec = do + describe "Special values" $ do + it "Infinity + 1 = Infinity" $ + op2 Op.add "Infinity" "1" `shouldBe` "Infinity" + it "NaN + 1 = NaN" $ + op2 Op.add "NaN" "1" `shouldBe` "NaN" + it "NaN + Infinity = NaN" $ + op2 Op.add "NaN" "Infinity" `shouldBe` "NaN" + it "1 - Infinity = -Infinity" $ + op2 Op.subtract "1" "Infinity" `shouldBe` "-Infinity" + it "-1 * Infinity = -Infinity" $ + op2 Op.multiply "-1" "Infinity" `shouldBe` "-Infinity" + it "-0 - 0 = -0" $ + op2 Op.subtract "-0" "0" `shouldBe` "-0" + it "-1 * 0 = -0" $ + op2 Op.multiply "-1" "0" `shouldBe` "-0" + it "1 / 0 = Infinity" $ + op2 Op.divide "1" "0" `shouldBe` "Infinity" + it "1 / -0 = -Infinity" $ + op2 Op.divide "1" "-0" `shouldBe` "-Infinity" + it "-1 / 0 = -Infinity" $ + op2 Op.divide "-1" "0" `shouldBe` "-Infinity" + + describe "add" $ do + it "add('12', '7.00') ==> '19.00'" $ + op2 Op.add "12" "7.00" `shouldBe` "19.00" + it "add('1E+2', '1E+4') ==> '1.01E+4'" $ + op2 Op.add "1E+2" "1E+4" `shouldBe` "1.01E+4" + + describe "subtract" $ do + it "subtract('1.3', '1.07') ==> '0.23'" $ + op2 Op.subtract "1.3" "1.07" `shouldBe` "0.23" + it "subtract('1.3', '1.30') ==> '0.00'" $ + op2 Op.subtract "1.3" "1.30" `shouldBe` "0.00" + it "subtract('1.3', '2.07') ==> '-0.77'" $ + op2 Op.subtract "1.3" "2.07" `shouldBe` "-0.77" + + describe "minus" $ do + it "minus('1.3') ==> '-1.3'" $ + op1 Op.minus "1.3" `shouldBe` "-1.3" + it "minus('-1.3') ==> '1.3'" $ + op1 Op.minus "-1.3" `shouldBe` "1.3" + + describe "plus" $ do + it "plus('1.3') ==> '1.3'" $ + op1 Op.plus "1.3" `shouldBe` "1.3" + it "plus('-1.3') ==> '-1.3'" $ + op1 Op.plus "-1.3" `shouldBe` "-1.3" + + describe "multiply" $ do + it "multiply('1.20', '3') ==> '3.60'" $ + op2 Op.multiply "1.20" "3" `shouldBe` "3.60" + it "multiply('7', '3') ==> '21'" $ + op2 Op.multiply "7" "3" `shouldBe` "21" + it "multiply('0.9', '0.8') ==> '0.72'" $ + op2 Op.multiply "0.9" "0.8" `shouldBe` "0.72" + it "multiply('0.9', '-0') ==> '-0.0'" $ + op2 Op.multiply "0.9" "-0" `shouldBe` "-0.0" + it "multiply('654321', '654321') ==> '4.28135971E+11'" $ + op2 Op.multiply "654321" "654321" `shouldBe` "4.28135971E+11" + + describe "exp" $ do + it "exp('-Infinity') ==> '0'" $ + op1 Op.exp "-Infinity" `shouldBe` "0" + it "exp('-1') ==> '0.367879441'" $ + op1 Op.exp "-1" `shouldBe` "0.367879441" + it "exp('0') ==> '1'" $ + op1 Op.exp "0" `shouldBe` "1" + it "exp('1') ==> '2.71828183'" $ + op1 Op.exp "1" `shouldBe` "2.71828183" + it "exp('0.693147181') ==> '2.00000000'" $ + op1 Op.exp "0.693147181" `shouldBe` "2.00000000" + it "exp('+Infinity') ==> 'Infinity'" $ + op1 Op.exp "+Infinity" `shouldBe` "Infinity" + + describe "fusedMultiplyAdd" $ do + it ("fused-multiply-add('3', '5', '7') ==> " ++ + "'22'") $ + op3 Op.fusedMultiplyAdd "3" "5" "7" `shouldBe` "22" + it ("fused-multiply-add('3', '-5', '7') ==> " ++ + "'-8'") $ + op3 Op.fusedMultiplyAdd "3" "-5" "7" `shouldBe` "-8" + it ("fused-multiply-add('888565290', '1557.96930', '-86087.7578') ==> " ++ + "'1.38435736E+12'") $ + op3 Op.fusedMultiplyAdd "888565290" "1557.96930" "-86087.7578" `shouldBe` + "1.38435736E+12" + + describe "ln" $ do + it "ln('0') ==> '-Infinity'" $ + op1 Op.ln "0" `shouldBe` "-Infinity" + it "ln('1.000') ==> '0'" $ + op1 Op.ln "1.000" `shouldBe` "0" + it "ln('2.71828183') ==> '1.00000000'" $ + op1 Op.ln "2.71828183" `shouldBe` "1.00000000" + it "ln('10') ==> '2.30258509'" $ + op1 Op.ln "10" `shouldBe` "2.30258509" + it "ln('+Infinity') ==> 'Infinity'" $ + op1 Op.ln "+Infinity" `shouldBe` "Infinity" + + describe "log10" $ do + it "log10('0') ==> '-Infinity'" $ + op1 Op.log10 "0" `shouldBe` "-Infinity" + it "log10('0.001') ==> '-3'" $ + op1 Op.log10 "0.001" `shouldBe` "-3" + it "log10('1.000') ==> '0'" $ + op1 Op.log10 "1.000" `shouldBe` "0" + it "log10('2') ==> '0.301029996'" $ + op1 Op.log10 "2" `shouldBe` "0.301029996" + it "log10('10') ==> '1'" $ + op1 Op.log10 "10" `shouldBe` "1" + it "log10('70') ==> '1.84509804'" $ + op1 Op.log10 "70" `shouldBe` "1.84509804" + it "log10('+Infinity') ==> 'Infinity'" $ + op1 Op.log10 "+Infinity" `shouldBe` "Infinity" + + describe "divide" $ do + it "divide('1', '3' ) ==> '0.333333333'" $ + op2 Op.divide "1" "3" `shouldBe` "0.333333333" + it "divide('2', '3' ) ==> '0.666666667'" $ + op2 Op.divide "2" "3" `shouldBe` "0.666666667" + it "divide('5', '2' ) ==> '2.5'" $ + op2 Op.divide "5" "2" `shouldBe` "2.5" + it "divide('1', '10' ) ==> '0.1'" $ + op2 Op.divide "1" "10" `shouldBe` "0.1" + it "divide('12', '12') ==> '1'" $ + op2 Op.divide "12" "12" `shouldBe` "1" + it "divide('8.00', '2') ==> '4.00'" $ + op2 Op.divide "8.00" "2" `shouldBe` "4.00" + it "divide('2.400', '2.0') ==> '1.20'" $ + op2 Op.divide "2.400" "2.0" `shouldBe` "1.20" + it "divide('1000', '100') ==> '10'" $ + op2 Op.divide "1000" "100" `shouldBe` "10" + it "divide('1000', '1') ==> '1000'" $ + op2 Op.divide "1000" "1" `shouldBe` "1000" + it "divide('2.40E+6', '2') ==> '1.20E+6'" $ + op2 Op.divide "2.40E+6" "2" `shouldBe` "1.20E+6" + + describe "abs" $ do + it "abs('2.1') ==> '2.1'" $ + op1 Op.abs "2.1" `shouldBe` "2.1" + it "abs('-100') ==> '100'" $ + op1 Op.abs "-100" `shouldBe` "100" + it "abs('101.5') ==> '101.5'" $ + op1 Op.abs "101.5" `shouldBe` "101.5" + it "abs('-101.5') ==> '101.5'" $ + op1 Op.abs "-101.5" `shouldBe` "101.5" + + describe "compare" $ do + it "compare('2.1', '3') ==> '-1'" $ + op2 compare' "2.1" "3" `shouldBe` "-1" + it "compare('2.1', '2.1') ==> '0'" $ + op2 compare' "2.1" "2.1" `shouldBe` "0" + it "compare('2.1', '2.10') ==> '0'" $ + op2 compare' "2.1" "2.10" `shouldBe` "0" + it "compare('3', '2.1') ==> '1'" $ + op2 compare' "3" "2.1" `shouldBe` "1" + it "compare('2.1', '-3') ==> '1'" $ + op2 compare' "2.1" "-3" `shouldBe` "1" + it "compare('-3', '2.1') ==> '-1'" $ + op2 compare' "-3" "2.1" `shouldBe` "-1" + + describe "max" $ do + it "max('3', '2') ==> '3'" $ + op2 Op.max "3" "2" `shouldBe` "3" + it "max('-10', '3') ==> '3'" $ + op2 Op.max "-10" "3" `shouldBe` "3" + it "max('1.0', '1') ==> '1'" $ + op2 Op.max "1.0" "1" `shouldBe` "1" + it "max('7', 'NaN') ==> '7'" $ + op2 Op.max "7" "NaN" `shouldBe` "7" + + describe "min" $ do + it "min('3', '2') ==> '2'" $ + op2 Op.min "3" "2" `shouldBe` "2" + it "min('-10', '3') ==> '-10'" $ + op2 Op.min "-10" "3" `shouldBe` "-10" + it "min('1.0', '1') ==> '1.0'" $ + op2 Op.min "1.0" "1" `shouldBe` "1.0" + it "min('7', 'NaN') ==> '7'" $ + op2 Op.min "7" "NaN" `shouldBe` "7" + + describe "power" $ do + it "power('2', '3') ==> '8'" $ + op2 Op.power "2" "3" `shouldBe` "8" + it "power('-2', '3') ==> '-8'" $ + op2 Op.power "-2" "3" `shouldBe` "-8" + it "power('2', '-3') ==> '0.125'" $ + op2 Op.power "2" "-3" `shouldBe` "0.125" + it "power('1.7', '8') ==> '69.7575744'" $ + op2 Op.power "1.7" "8" `shouldBe` "69.7575744" + it "power('10', '0.301029996') ==> '2.00000000'" $ + op2 Op.power "10" "0.301029996" `shouldBe` "2.00000000" + it "power('Infinity', '-1') ==> '0'" $ + op2 Op.power "Infinity" "-1" `shouldBe` "0" + it "power('Infinity', '0') ==> '1'" $ + op2 Op.power "Infinity" "0" `shouldBe` "1" + it "power('Infinity', '1') ==> 'Infinity'" $ + op2 Op.power "Infinity" "1" `shouldBe` "Infinity" + it "power('-Infinity', '-1') ==> '-0'" $ + op2 Op.power "-Infinity" "-1" `shouldBe` "-0" + it "power('-Infinity', '0') ==> '1'" $ + op2 Op.power "-Infinity" "0" `shouldBe` "1" + it "power('-Infinity', '1') ==> '-Infinity'" $ + op2 Op.power "-Infinity" "1" `shouldBe` "-Infinity" + it "power('-Infinity', '2') ==> 'Infinity'" $ + op2 Op.power "-Infinity" "2" `shouldBe` "Infinity" + it "power('0', '0') ==> 'NaN'" $ + op2 Op.power "0" "0" `shouldBe` "NaN" + + describe "quantize" $ do + it "quantize('2.17', '0.001') ==> '2.170'" $ + op2 Op.quantize "2.17" "0.001" `shouldBe` "2.170" + it "quantize('2.17', '0.01') ==> '2.17'" $ + op2 Op.quantize "2.17" "0.01" `shouldBe` "2.17" + it "quantize('2.17', '0.1') ==> '2.2'" $ + op2 Op.quantize "2.17" "0.1" `shouldBe` "2.2" + it "quantize('2.17', '1e+0') ==> '2'" $ + op2 Op.quantize "2.17" "1e+0" `shouldBe` "2" + it "quantize('2.17', '1e+1') ==> '0E+1'" $ + op2 Op.quantize "2.17" "1e+1" `shouldBe` "0E+1" + it "quantize('-Inf' 'Infinity') ==> '-Infinity'" $ + op2 Op.quantize "-Inf" "Infinity" `shouldBe` "-Infinity" + it "quantize('2', 'Infinity') ==> 'NaN'" $ + op2 Op.quantize "2" "Infinity" `shouldBe` "NaN" + it "quantize('-0.1', '1' ) ==> '-0'" $ + op2 Op.quantize "-0.1" "1" `shouldBe` "-0" + it "quantize('-0', '1e+5') ==> '-0E+5'" $ + op2 Op.quantize "-0" "1e+5" `shouldBe` "-0E+5" + it "quantize('+35236450.6', '1e-2') ==> 'NaN'" $ + op2 Op.quantize "+35236450.6" "1e-2" `shouldBe` "NaN" + it "quantize('-35236450.6', '1e-2') ==> 'NaN'" $ + op2 Op.quantize "-35236450.6" "1e-2" `shouldBe` "NaN" + it "quantize('217', '1e-1') ==> '217.0'" $ + op2 Op.quantize "217" "1e-1" `shouldBe` "217.0" + it "quantize('217', '1e+0') ==> '217'" $ + op2 Op.quantize "217" "1e+0" `shouldBe` "217" + it "quantize('217', '1e+1') ==> '2.2E+2'" $ + op2 Op.quantize "217" "1e+1" `shouldBe` "2.2E+2" + it "quantize('217', '1e+2') ==> '2E+2'" $ + op2 Op.quantize "217" "1e+2" `shouldBe` "2E+2" + + describe "reduce" $ do + it "reduce('2.1') ==> '2.1'" $ + op1 Op.reduce "2.1" `shouldBe` "2.1" + it "reduce('-2.0') ==> '-2'" $ + op1 Op.reduce "-2.0" `shouldBe` "-2" + it "reduce('1.200') ==> '1.2'" $ + op1 Op.reduce "1.200" `shouldBe` "1.2" + it "reduce('-120') ==> '-1.2E+2'" $ + op1 Op.reduce "-120" `shouldBe` "-1.2E+2" + it "reduce('120.00') ==> '1.2E+2'" $ + op1 Op.reduce "120.00" `shouldBe` "1.2E+2" + it "reduce('0.00') ==> '0'" $ + op1 Op.reduce "0.00" `shouldBe` "0" + + describe "roundToIntegralExact" $ do + it "round-to-integral-exact('2.1') ==> '2'" $ + op1 Op.roundToIntegralExact "2.1" `shouldBe` "2" + it "round-to-integral-exact('100') ==> '100'" $ + op1 Op.roundToIntegralExact "100" `shouldBe` "100" + it "round-to-integral-exact('100.0') ==> '100'" $ + op1 Op.roundToIntegralExact "100.0" `shouldBe` "100" + it "round-to-integral-exact('101.5') ==> '102'" $ + op1 Op.roundToIntegralExact "101.5" `shouldBe` "102" + it "round-to-integral-exact('-101.5') ==> '-102'" $ + op1 Op.roundToIntegralExact "-101.5" `shouldBe` "-102" + it "round-to-integral-exact('10E+5') ==> '1.0E+6'" $ + op1 Op.roundToIntegralExact "10E+5" `shouldBe` "1.0E+6" + it "round-to-integral-exact('7.89E+77') ==> '7.89E+77'" $ + op1 Op.roundToIntegralExact "7.89E+77" `shouldBe` "7.89E+77" + it "round-to-integral-exact('-Inf') ==> '-Infinity'" $ + op1 Op.roundToIntegralExact "-Inf" `shouldBe` "-Infinity" + + describe "squareRoot" $ do + it "square-root('0') ==> '0'" $ + op1 Op.squareRoot "0" `shouldBe` "0" + it "square-root('-0') ==> '-0'" $ + op1 Op.squareRoot "-0" `shouldBe` "-0" + -- The following example is a corrected version of that found in the + -- specification; confirmed with Mike Cowlishaw on 2016-08-02. + it "square-root('0.39') ==> '0.624499800'" $ + op1 Op.squareRoot "0.39" `shouldBe` "0.624499800" + it "square-root('100') ==> '10'" $ + op1 Op.squareRoot "100" `shouldBe` "10" + it "square-root('1') ==> '1'" $ + op1 Op.squareRoot "1" `shouldBe` "1" + it "square-root('1.0') ==> '1.0'" $ + op1 Op.squareRoot "1.0" `shouldBe` "1.0" + it "square-root('1.00') ==> '1.0'" $ + op1 Op.squareRoot "1.00" `shouldBe` "1.0" + it "square-root('7') ==> '2.64575131'" $ + op1 Op.squareRoot "7" `shouldBe` "2.64575131" + it "square-root('10') ==> '3.16227766'" $ + op1 Op.squareRoot "10" `shouldBe` "3.16227766" + + describe "and" $ do + it "and('0', '0') ==> '0'" $ + op2 Op.and "0" "0" `shouldBe` "0" + it "and('0', '1') ==> '0'" $ + op2 Op.and "0" "1" `shouldBe` "0" + it "and('1', '0') ==> '0'" $ + op2 Op.and "1" "0" `shouldBe` "0" + it "and('1', '1') ==> '1'" $ + op2 Op.and "1" "1" `shouldBe` "1" + it "and('1100', '1010') ==> '1000'" $ + op2 Op.and "1100" "1010" `shouldBe` "1000" + it "and('1111', '10') ==> '10'" $ + op2 Op.and "1111" "10" `shouldBe` "10" + + describe "or" $ do + it "or('0', '0') ==> '0'" $ + op2 Op.or "0" "0" `shouldBe` "0" + it "or('0', '1') ==> '1'" $ + op2 Op.or "0" "1" `shouldBe` "1" + it "or('1', '0') ==> '1'" $ + op2 Op.or "1" "0" `shouldBe` "1" + it "or('1', '1') ==> '1'" $ + op2 Op.or "1" "1" `shouldBe` "1" + it "or('1100', '1010') ==> '1110'" $ + op2 Op.or "1100" "1010" `shouldBe` "1110" + it "or('1110', '10') ==> '1110'" $ + op2 Op.or "1110" "10" `shouldBe` "1110" + + describe "xor" $ do + it "xor('0', '0') ==> '0'" $ + op2 Op.xor "0" "0" `shouldBe` "0" + it "xor('0', '1') ==> '1'" $ + op2 Op.xor "0" "1" `shouldBe` "1" + it "xor('1', '0') ==> '1'" $ + op2 Op.xor "1" "0" `shouldBe` "1" + it "xor('1', '1') ==> '0'" $ + op2 Op.xor "1" "1" `shouldBe` "0" + it "xor('1100', '1010') ==> '110'" $ + op2 Op.xor "1100" "1010" `shouldBe` "110" + it "xor('1111', '10') ==> '1101'" $ + op2 Op.xor "1111" "10" `shouldBe` "1101" + + describe "invert" $ do + it "invert('0') ==> '111111111'" $ + op1 Op.invert "0" `shouldBe` "111111111" + it "invert('1') ==> '111111110'" $ + op1 Op.invert "1" `shouldBe` "111111110" + it "invert('111111111') ==> '0'" $ + op1 Op.invert "111111111" `shouldBe` "0" + it "invert('101010101') ==> '10101010'" $ + op1 Op.invert "101010101" `shouldBe` "10101010" + + describe "canonical" $ + it "canonical('2.50') ==> '2.50'" $ + op1 Op.canonical "2.50" `shouldBe` "2.50" + + describe "class'" $ do + it "class('Infinity') ==> \"+Infinity\"" $ + op1 Op.class' "Infinity" `shouldBe` "+Infinity" + it "class('1E-10') ==> \"+Normal\"" $ + op1 Op.class' "1E-10" `shouldBe` "+Normal" + it "class('2.50') ==> \"+Normal\"" $ + op1 Op.class' "2.50" `shouldBe` "+Normal" + it "class('0.1E-999') ==> \"+Subnormal\"" $ + op1 Op.class' "0.1E-999" `shouldBe` "+Subnormal" + it "class('0') ==> \"+Zero\"" $ + op1 Op.class' "0" `shouldBe` "+Zero" + it "class('-0') ==> \"-Zero\"" $ + op1 Op.class' "-0" `shouldBe` "-Zero" + it "class('-0.1E-999') ==> \"-Subnormal\"" $ + op1 Op.class' "-0.1E-999" `shouldBe` "-Subnormal" + it "class('-1E-10') ==> \"-Normal\"" $ + op1 Op.class' "-1E-10" `shouldBe` "-Normal" + it "class('-2.50') ==> \"-Normal\"" $ + op1 Op.class' "-2.50" `shouldBe` "-Normal" + it "class('-Infinity') ==> \"-Infinity\"" $ + op1 Op.class' "-Infinity" `shouldBe` "-Infinity" + it "class('NaN') ==> \"NaN\"" $ + op1 Op.class' "NaN" `shouldBe` "NaN" + it "class('-NaN') ==> \"NaN\"" $ + op1 Op.class' "-NaN" `shouldBe` "NaN" + it "class('sNaN') ==> \"sNaN\"" $ + op1 Op.class' "sNaN" `shouldBe` "sNaN" + + describe "compareTotal" $ do + it "compare-total('12.73', '127.9') ==> '-1'" $ + op2 compareTotal' "12.73" "127.9" `shouldBe` "-1" + it "compare-total('-127', '12') ==> '-1'" $ + op2 compareTotal' "-127" "12" `shouldBe` "-1" + it "compare-total('12.30', '12.3') ==> '-1'" $ + op2 compareTotal' "12.30" "12.3" `shouldBe` "-1" + it "compare-total('12.30', '12.30') ==> '0'" $ + op2 compareTotal' "12.30" "12.30" `shouldBe` "0" + it "compare-total('12.3', '12.300') ==> '1'" $ + op2 compareTotal' "12.3" "12.300" `shouldBe` "1" + it "compare-total('12.3', 'NaN') ==> '-1'" $ + op2 compareTotal' "12.3" "NaN" `shouldBe` "-1" + + describe "copy" $ do + it "copy('2.1') ==> '2.1'" $ + op1 Op.copy "2.1" `shouldBe` "2.1" + it "copy('-1.00') ==> '-1.00'" $ + op1 Op.copy "-1.00" `shouldBe` "-1.00" + + describe "copyAbs" $ do + it "copy-abs('2.1') ==> '2.1'" $ + op1 Op.copyAbs "2.1" `shouldBe` "2.1" + it "copy-abs('-100') ==> '100'" $ + op1 Op.copyAbs "-100" `shouldBe` "100" + + describe "copyNegate" $ do + it "copy-negate('101.5') ==> '-101.5'" $ + op1 Op.copyNegate "101.5" `shouldBe` "-101.5" + it "copy-negate('-101.5') ==> '101.5'" $ + op1 Op.copyNegate "-101.5" `shouldBe` "101.5" + + describe "copySign" $ do + it "copy-sign( '1.50', '7.33') ==> '1.50'" $ + op2 Op.copySign "1.50" "7.33" `shouldBe` "1.50" + it "copy-sign('-1.50', '7.33') ==> '1.50'" $ + op2 Op.copySign "-1.50" "7.33" `shouldBe` "1.50" + it "copy-sign( '1.50', '-7.33') ==> '-1.50'" $ + op2 Op.copySign "1.50" "-7.33" `shouldBe` "-1.50" + it "copy-sign('-1.50', '-7.33') ==> '-1.50'" $ + op2 Op.copySign "-1.50" "-7.33" `shouldBe` "-1.50" + + describe "isCanonical" $ + it "is-canonical('2.50') ==> '1'" $ + pred1 Op.isCanonical "2.50" `shouldBe` "1" + + describe "isFinite" $ do + it "is-finite('2.50') ==> '1'" $ + pred1 Op.isFinite "2.50" `shouldBe` "1" + it "is-finite('-0.3') ==> '1'" $ + pred1 Op.isFinite "-0.3" `shouldBe` "1" + it "is-finite('0') ==> '1'" $ + pred1 Op.isFinite "0" `shouldBe` "1" + it "is-finite('Inf') ==> '0'" $ + pred1 Op.isFinite "Inf" `shouldBe` "0" + it "is-finite('NaN') ==> '0'" $ + pred1 Op.isFinite "NaN" `shouldBe` "0" + + describe "isInfinite" $ do + it "is-infinite('2.50') ==> '0'" $ + pred1 Op.isInfinite "2.50" `shouldBe` "0" + it "is-infinite('-Inf') ==> '1'" $ + pred1 Op.isInfinite "-Inf" `shouldBe` "1" + it "is-infinite('NaN') ==> '0'" $ + pred1 Op.isInfinite "NaN" `shouldBe` "0" + + describe "isNaN" $ do + it "is-NaN('2.50') ==> '0'" $ + pred1 Op.isNaN "2.50" `shouldBe` "0" + it "is-NaN('NaN') ==> '1'" $ + pred1 Op.isNaN "NaN" `shouldBe` "1" + it "is-NaN('-sNaN') ==> '1'" $ + pred1 Op.isNaN "-sNaN" `shouldBe` "1" + + describe "isNormal" $ do + it "is-normal('2.50') ==> '1'" $ + pred1 Op.isNormal "2.50" `shouldBe` "1" + it "is-normal('0.1E-999') ==> '0'" $ + pred1 Op.isNormal "0.1E-999" `shouldBe` "0" + it "is-normal('0.00') ==> '0'" $ + pred1 Op.isNormal "0.00" `shouldBe` "0" + it "is-normal('-Inf') ==> '0'" $ + pred1 Op.isNormal "-Inf" `shouldBe` "0" + it "is-normal('NaN') ==> '0'" $ + pred1 Op.isNormal "NaN" `shouldBe` "0" + + describe "isQNaN" $ do + it "is-qNaN('2.50') ==> '0'" $ + pred1 Op.isQNaN "2.50" `shouldBe` "0" + it "is-qNaN('NaN') ==> '1'" $ + pred1 Op.isQNaN "NaN" `shouldBe` "1" + it "is-qNaN('sNaN') ==> '0'" $ + pred1 Op.isQNaN "sNaN" `shouldBe` "0" + + describe "isSigned" $ do + it "is-signed('2.50') ==> '0'" $ + pred1 Op.isSigned "2.50" `shouldBe` "0" + it "is-signed('-12') ==> '1'" $ + pred1 Op.isSigned "-12" `shouldBe` "1" + it "is-signed('-0') ==> '1'" $ + pred1 Op.isSigned "-0" `shouldBe` "1" + + describe "isSNaN" $ do + it "is-sNaN('2.50') ==> '0'" $ + pred1 Op.isSNaN "2.50" `shouldBe` "0" + it "is-sNaN('NaN') ==> '0'" $ + pred1 Op.isSNaN "NaN" `shouldBe` "0" + it "is-sNaN('sNaN') ==> '1'" $ + pred1 Op.isSNaN "sNaN" `shouldBe` "1" + + describe "isSubnormal" $ do + it "is-subnormal('2.50') ==> '0'" $ + pred1 Op.isSubnormal "2.50" `shouldBe` "0" + it "is-subnormal('0.1E-999') ==> '1'" $ + pred1 Op.isSubnormal "0.1E-999" `shouldBe` "1" + it "is-subnormal('0.00') ==> '0'" $ + pred1 Op.isSubnormal "0.00" `shouldBe` "0" + it "is-subnormal('-Inf') ==> '0'" $ + pred1 Op.isSubnormal "-Inf" `shouldBe` "0" + it "is-subnormal('NaN') ==> '0'" $ + pred1 Op.isSubnormal "NaN" `shouldBe` "0" + + describe "isZero" $ do + it "is-zero('0') ==> '1'" $ + pred1 Op.isZero "0" `shouldBe` "1" + it "is-zero('2.50') ==> '0'" $ + pred1 Op.isZero "2.50" `shouldBe` "0" + it "is-zero('-0E+2') ==> '1'" $ + pred1 Op.isZero "-0E+2" `shouldBe` "1" + + describe "logb" $ do + it "logb('250') ==> '2'" $ + op1 Op.logb "250" `shouldBe` "2" + it "logb('2.50') ==> '0'" $ + op1 Op.logb "2.50" `shouldBe` "0" + it "logb('0.03') ==> '-2'" $ + op1 Op.logb "0.03" `shouldBe` "-2" + it "logb('0') ==> '-Infinity'" $ + op1 Op.logb "0" `shouldBe` "-Infinity" + + describe "scaleb" $ do + it "scaleb('7.50', '-2') ==> '0.0750'" $ + op2 Op.scaleb "7.50" "-2" `shouldBe` "0.0750" + it "scaleb('7.50', '0') ==> '7.50'" $ + op2 Op.scaleb "7.50" "0" `shouldBe` "7.50" + it "scaleb('7.50', '3') ==> '7.50E+3'" $ + op2 Op.scaleb "7.50" "3" `shouldBe` "7.50E+3" + + describe "radix" $ + it "radix() ==> '10'" $ + op0 Op.radix `shouldBe` "10" + + describe "sameQuantum" $ do + it "samequantum('2.17', '0.001') ==> '0'" $ + pred2 Op.sameQuantum "2.17" "0.001" `shouldBe` "0" + it "samequantum('2.17', '0.01') ==> '1'" $ + pred2 Op.sameQuantum "2.17" "0.01" `shouldBe` "1" + it "samequantum('2.17', '0.1') ==> '0'" $ + pred2 Op.sameQuantum "2.17" "0.1" `shouldBe` "0" + it "samequantum('2.17', '1') ==> '0'" $ + pred2 Op.sameQuantum "2.17" "1" `shouldBe` "0" + it "samequantum('Inf', '-Inf') ==> '1'" $ + pred2 Op.sameQuantum "Inf" "-Inf" `shouldBe` "1" + it "samequantum('NaN', 'NaN') ==> '1'" $ + pred2 Op.sameQuantum "NaN" "NaN" `shouldBe` "1" + + describe "shift" $ do + it "shift('34', '8') ==> '400000000'" $ + op2 Op.shift "34" "8" `shouldBe` "400000000" + it "shift('12', '9') ==> '0'" $ + op2 Op.shift "12" "9" `shouldBe` "0" + it "shift('123456789', '-2') ==> '1234567'" $ + op2 Op.shift "123456789" "-2" `shouldBe` "1234567" + it "shift('123456789', '0') ==> '123456789'" $ + op2 Op.shift "123456789" "0" `shouldBe` "123456789" + it "shift('123456789', '+2') ==> '345678900'" $ + op2 Op.shift "123456789" "+2" `shouldBe` "345678900" + + describe "rotate" $ do + it "rotate('34', '8') ==> '400000003'" $ + op2 Op.rotate "34" "8" `shouldBe` "400000003" + it "rotate('12', '9') ==> '12'" $ + op2 Op.rotate "12" "9" `shouldBe` "12" + it "rotate('123456789', '-2') ==> '891234567'" $ + op2 Op.rotate "123456789" "-2" `shouldBe` "891234567" + it "rotate('123456789', '0') ==> '123456789'" $ + op2 Op.rotate "123456789" "0" `shouldBe` "123456789" + it "rotate('123456789', '+2') ==> '345678912'" $ + op2 Op.rotate "123456789" "+2" `shouldBe` "345678912" + +exceptionError :: Exception p r -> a +exceptionError = error . show . exceptionSignal + +type BasicArith = Arith P9 RoundHalfUp + +pred1 :: (BasicDecimal -> BasicArith Bool) -> String -> String +pred1 op x = either exceptionError (show . fromBool) $ + evalArith arith newContext + where arith = op (read x) + +pred2 :: (BasicDecimal -> BasicDecimal -> BasicArith Bool) -> String -> String + -> String +pred2 op x y = either exceptionError (show . fromBool) $ + evalArith arith newContext + where arith = op (read x) (read y) + +op0 :: Show a => BasicArith a -> String +op0 op = either exceptionError show $ evalArith op newContext + +op1 :: Show a => (BasicDecimal -> BasicArith a) -> String -> String +op1 op x = either exceptionError show $ evalArith arith newContext + where arith = op (read x) + +op2 :: Show a => (BasicDecimal -> BasicDecimal -> BasicArith a) -> String + -> String -> String +op2 op x y = either exceptionError show $ evalArith arith newContext + where arith = read x `op` read y + +op3 :: Show a => (BasicDecimal -> BasicDecimal -> BasicDecimal -> BasicArith a) + -> String -> String -> String -> String +op3 op x y z = either exceptionError show $ evalArith arith newContext + where arith = op (read x) (read y) (read z) + +compare' :: Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +compare' x y = either id fromOrdering <$> Op.compare x y + +compareTotal' :: Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +compareTotal' x y = fromOrdering <$> Op.compareTotal x y diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -0,0 +1,1 @@ +{-# OPTIONS_GHC -F -pgmF hspec-discover #-}