diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -11,4 +11,4 @@ [General Decimal Arithmetic]: http://speleotrove.com/decimal/ While usable, the implementation is currently in its infancy. Additional -operations as well as an API for manipulating context flags are planned. +operations and possible API changes are planned. diff --git a/TODO b/TODO --- a/TODO +++ b/TODO @@ -1,7 +1,5 @@ -*- Outline -*- * To Do -** Test suite -** instance Floating (Number p r) -** instance RealFloat (Number p r) -** instance PrintfArg (Number p r) +** 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.1.0.1 +version: 0.2.0.0 synopsis: An implementation of Mike Cowlishaw's General Decimal Arithmetic Specification @@ -39,22 +39,27 @@ exposed-modules: Numeric.Decimal Numeric.Decimal.Conversion + Numeric.Decimal.Arithmetic Numeric.Decimal.Operation other-modules: Numeric.Decimal.Number Numeric.Decimal.Precision Numeric.Decimal.Rounding build-depends: base >= 4.7 && < 5 + , mtl default-language: Haskell2010 default-extensions: Trustworthy - other-extensions: RoleAnnotations + other-extensions: FlexibleInstances + MultiParamTypeClasses + RoleAnnotations -test-suite decimal-arithmetic-test +test-suite doctests type: exitcode-stdio-1.0 hs-source-dirs: test - main-is: Spec.hs + main-is: doctests.hs build-depends: base , decimal-arithmetic + , doctest >= 0.8 , QuickCheck ghc-options: -threaded -rtsopts -with-rtsopts=-N default-language: Haskell2010 diff --git a/src/Numeric/Decimal.hs b/src/Numeric/Decimal.hs --- a/src/Numeric/Decimal.hs +++ b/src/Numeric/Decimal.hs @@ -6,7 +6,7 @@ Maintainer : rob@mars.org Stability : experimental -This module provides a general-purpose 'Number' type supporting decimal +This module provides a general-purpose number type supporting decimal arithmetic for both limited precision floating-point (IEEE 754-2008) and for arbitrary precision floating-point (following the same principles as IEEE 754 and IEEE 854-1987) as described in the @@ -14,13 +14,28 @@ by Mike Cowlishaw. In addition to floating-point arithmetic, integer and unrounded floating-point arithmetic are included as subsets. -Unlike the binary floating-point types 'Float' and 'Double', the 'Number' type -can represent and perform arithmetic with decimal numbers exactly. -Internally, a 'Number' is represented with an integral coefficient and base-10 -exponent. +Unlike the binary floating-point types 'Float' and 'Double', decimal number +types can perform decimal arithmetic exactly. Internally, decimal numbers are +represented with an integral coefficient and base-10 exponent. -The 'Number' type supports lossless conversion to and from a string -representation via the 'Show' and 'Read' instances. Note that there may be +@ +@>>>@ __29.99 + 4.71 :: Double__ +34.699999999999996 +@ + +>>> 29.99 + 4.71 :: BasicDecimal +34.70 + +@ +@>>>@ __0.1 + 0.2 == (0.3 :: Double)__ +False +@ + +>>> 0.1 + 0.2 == (0.3 :: BasicDecimal) +True + +Decimal numbers support lossless conversion to and from a string +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. -} @@ -29,7 +44,7 @@ -- $usage -- * Arbitrary-precision decimal numbers - Number + Decimal , BasicDecimal , ExtendedDecimal , GeneralDecimal @@ -51,61 +66,61 @@ -- * Functions , cast + , toBool ) where import Numeric.Decimal.Number import Numeric.Decimal.Precision import Numeric.Decimal.Rounding --- | A basic decimal floating point number with 9 digits of precision, rounding half up -type BasicDecimal = Number P9 RoundHalfUp +-- | 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 = Number p RoundHalfEven +-- | 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 -basicDefaultContext :: TrapHandler P9 RoundHalfUp -> Context P9 RoundHalfUp -basicDefaultContext handler = defaultContext { trapHandler = trap } - where trap Inexact = id - trap Rounded = id - trap Subnormal = id - trap sig = handler sig +{- $usage -extendedDefaultContext :: Context p RoundHalfEven -extendedDefaultContext = defaultContext +You should choose a decimal number type with appropriate precision and +rounding to use in your application. There are several options: --- $usage --- --- It is recommended to create an alias for the type of numbers you wish to --- support in your application. For example: --- --- > type Decimal = BasicDecimal --- --- This is a basic number type with 9 decimal digits of precision that rounds --- half up. --- --- > type Decimal = ExtendedDecimal P15 --- --- This is a number type with 15 decimal digits of precision that rounds half --- even. There are a range of ready-made precisions available, including 'P1' --- through 'P50' on up to 'P2000'. Alternatively, an arbitrary precision can --- be constructed through type application of 'PPlus1' and/or 'PTimes2' to any --- existing precision. --- --- > type Decimal = GeneralDecimal --- --- This is a number type with infinite precision. Note that not all operations --- support numbers with infinite precision. --- --- > type Decimal = Number P34 RoundDown --- --- This is a custom number type with 34 decimal digits of precision that --- rounds down (truncates). Several 'Rounding' algorithms are available to --- choose from. --- --- A decimal number type may be used in a @default@ declaration, possibly --- replacing 'Double' or 'Integer'. For example: --- --- > default (Integer, Decimal) +* 'BasicDecimal' is a number type with 9 decimal digits of precision that +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. + +* 'GeneralDecimal' is a number type with infinite precision. Note that not all +operations support numbers with infinite precision. + +* 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: + +> default (Integer, BasicDecimal) + +== 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 +"Numeric.Decimal.Operation". +-} diff --git a/src/Numeric/Decimal/Arithmetic.hs b/src/Numeric/Decimal/Arithmetic.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Decimal/Arithmetic.hs @@ -0,0 +1,264 @@ + +{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-} + +-- | It is not usually necessary to import this module unless you want to use +-- the arithmetic operations from "Numeric.Decimal.Operation" or you need +-- precise control over the handling of exceptional conditions in an +-- arithmetic computation. + +module Numeric.Decimal.Arithmetic + ( -- * Decimal arithmetic + -- $decimal-arithmetic + + -- ** Context + Context + , newContext + , flags + + -- *** Default contexts + -- $default-contexts + , basicDefaultContext + , extendedDefaultContext + + -- ** The Arith monad + , Arith + , runArith + , evalArith + + -- * Exceptional conditions + -- $exceptional-conditions + , Exception + , exceptionSignal + , exceptionResult + + -- ** Signals + , Signal(..) + , Signals + , signal + , signals + , signalMember + + , raiseSignal + , clearFlags + + -- ** Traps + , TrapHandler + , trap + ) where + +import Control.Monad.Except (MonadError(throwError, catchError), + ExceptT, runExceptT) +import Control.Monad.State (MonadState(get, put), modify, gets, + State, runState, evalState) +import Data.Bits (bit, complement, testBit, (.&.), (.|.)) +import Data.Monoid ((<>)) + +import Numeric.Decimal.Number +import Numeric.Decimal.Precision +import Numeric.Decimal.Rounding + +-- $decimal-arithmetic +-- +-- Decimal arithmetic is performed within a context that maintains state to +-- handle exceptional conditions such as underflow, rounding, or division by +-- zero (cf. 'Signal'). The 'Arith' monad provides a means to evaluate an +-- 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@ +data Context p r = + Context { flags :: Signals + -- ^ The current signal flags of the context + , trapHandler :: TrapHandler p r + -- ^ The trap handler function for the context + } + +-- | Return a new context with all signal flags cleared and all traps disabled. +newContext :: Context p r +newContext = Context { flags = mempty + , trapHandler = return . exceptionResult + } + +-- $default-contexts +-- +-- The /General Decimal Arithmetic/ specification defines optional default +-- contexts, which define suitable settings for basic arithmetic and for the +-- extended arithmetic defined by IEEE 854 and IEEE 754. + +-- | Return a new context with all signal flags cleared, all traps enabled +-- except for 'Inexact', 'Rounded', and 'Subnormal', using a precision of 9 +-- significant decimal digits, and rounding half up. Trapped signals simply +-- call 'throwError' with the corresponding 'Exception', and can be caught +-- using 'catchError'. +basicDefaultContext :: Context P9 RoundHalfUp +basicDefaultContext = newContext { trapHandler = handler } + where handler e + | exceptionSignal e `notElem` disabled = throwError e + | otherwise = trapHandler newContext e + 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). +extendedDefaultContext :: Context p RoundHalfEven +extendedDefaultContext = newContext + +-- | A representation of an exceptional condition +data Exception p r = + Exception { exceptionSignal :: Signal + -- ^ The signal raised by the exceptional + -- condition + , exceptionResult :: Decimal p r + -- ^ The defined result for the exceptional + -- condition + } + deriving Show + +-- | A decimal arithmetic monad parameterized by the precision @p@ and +-- rounding algorithm @r@ +newtype Arith p r a = Arith (ExceptT (Exception p r) + (State (Context p r)) a) + +instance Functor (Arith p r) where + fmap f (Arith s) = Arith (fmap f s) + +instance Applicative (Arith p r) where + pure = Arith . pure + Arith f <*> Arith e = Arith (f <*> e) + +instance Monad (Arith p r) where + Arith e >>= f = Arith (e >>= g) + where g x = let Arith t = f x in t + +instance MonadError (Exception p r) (Arith p r) where + throwError = Arith . throwError + catchError (Arith e) f = Arith (catchError e g) + where g x = let Arith t = f x in t + +instance MonadState (Context p r) (Arith p r) where + get = Arith get + put = Arith . put + +-- | Evaluate an arithmetic computation in the given context and return the +-- final value (or exception) and resulting context. +runArith :: Arith p r a -> Context p r + -> (Either (Exception p r) a, Context p r) +runArith (Arith e) = runState (runExceptT e) + +-- | Evaluate an arithmetic computation in the given context and return the +-- final value or exception, discarding the resulting context. +evalArith :: Arith p r a -> Context p r -> Either (Exception p r) a +evalArith (Arith e) = evalState (runExceptT e) + +-- $exceptional-conditions +-- +-- Exceptional conditions are grouped into signals, which can be controlled +-- individually. A 'Context' contains a flag and a trap enabler (i.e. enabled +-- or disabled) for each 'Signal'. + +data Signal + = Clamped + -- ^ Raised when the exponent of a result has been altered or constrained + -- in order to fit the constraints of a specific concrete representation + | DivisionByZero + -- ^ Raised when a non-zero dividend is divided by zero + | Inexact + -- ^ Raised when a result is not exact (one or more non-zero coefficient + -- digits were discarded during rounding) + | InvalidOperation + -- ^ Raised when a result would be undefined or impossible + | Overflow + -- ^ Raised when the exponent of a result is too large to be represented + | Rounded + -- ^ Raised when a result has been rounded (that is, some zero or non-zero + -- coefficient digits were discarded) + | Subnormal + -- ^ Raised when a result is subnormal (its adjusted exponent is less than + -- E/min/), before any rounding + | Underflow + -- ^ Raised when a result is both subnormal and inexact + deriving (Eq, Enum, Bounded, Show) + +-- | A group of signals can be manipulated as a set. +newtype Signals = Signals Int + deriving Eq + +instance Show Signals where + showsPrec d sigs = showParen (d > 10) $ + showString "signals " . showsPrec 11 (signalList sigs) + +instance Monoid Signals where + mempty = Signals 0 + Signals x `mappend` Signals y = Signals (x .|. y) + +-- | Create a set of signals from a singleton. +signal :: Signal -> Signals +signal = Signals . bit . fromEnum + +-- | Create a set of signals from a list. +signals :: [Signal] -> Signals +signals = mconcat . map signal + +-- | Enumerate the given set of signals. +signalList :: Signals -> [Signal] +signalList sigs = filter (`signalMember` sigs) [minBound..maxBound] + +-- | Remove the first set of signals from the second. +unsignal :: Signals -> Signals -> Signals +unsignal (Signals u) (Signals ss) = Signals (ss .&. complement u) + +-- | Determine whether a signal is a member of a set. +signalMember :: Signal -> Signals -> Bool +signalMember sig (Signals ss) = testBit ss (fromEnum sig) + +-- | Set the given signal flag in the context of the current arithmetic +-- computation, and call the trap handler if the trap for this signal is +-- currently enabled. +raiseSignal :: Signal -> Decimal p r -> Arith p r (Decimal p r) +raiseSignal sig n = do + ctx <- get + let ctx' = ctx { flags = flags ctx <> signal sig } + put ctx' + trapHandler ctx' (Exception sig n) + +-- | Clear the given signal flags from the context of the current arithmetic +-- computation. +clearFlags :: Signals -> Arith p r () +clearFlags sigs = modify $ \ctx -> ctx { flags = unsignal sigs (flags ctx) } + +-- | A trap handler function may return a substitute result for the operation +-- that caused the exceptional condition, or it may call 'throwError' to abort +-- the arithmetic computation (or pass control to an enclosing 'catchError' +-- handler). +type TrapHandler p r = Exception p r -> Arith p r (Decimal p r) + +-- | Evaluate an arithmetic computation within a modified context that enables +-- the given signals to be trapped by the given handler. The previous trap +-- handler (and enabler state) will be restored during any trap, as well as +-- upon completion. Any existing trap handlers for signals not mentioned +-- remain in effect. +trap :: Signals -> TrapHandler p r -> Arith p r a -> Arith p r a +trap sigs handler arith = do + origHandler <- gets trapHandler + + let newHandler e = wrapHandler origHandler $ + if exceptionSignal e `signalMember` sigs + then handler e + else origHandler e + + wrapHandler newHandler arith `catchError` \e -> do + setHandler origHandler + throwError e + + where wrapHandler :: TrapHandler p r -> Arith p r a -> Arith p r a + wrapHandler handler arith = do + prevHandler <- gets trapHandler + setHandler handler + r <- arith + setHandler prevHandler + return r + + setHandler :: TrapHandler p r -> Arith p r () + setHandler handler = modify $ \ctx -> ctx { trapHandler = handler } diff --git a/src/Numeric/Decimal/Arithmetic.hs-boot b/src/Numeric/Decimal/Arithmetic.hs-boot new file mode 100644 --- /dev/null +++ b/src/Numeric/Decimal/Arithmetic.hs-boot @@ -0,0 +1,43 @@ +-- -*- Haskell -*- + +{-# LANGUAGE RoleAnnotations #-} + +module Numeric.Decimal.Arithmetic + ( Arith + , newContext + , evalArith + , Signal(..) + , raiseSignal + , exceptionResult + ) where + +import Control.Monad.Except (ExceptT) +import Control.Monad.State (State) + +import {-# SOURCE #-} Numeric.Decimal.Number + +--instance Precision p => Precision (Arith p r a) + +data Context p r +newContext :: Context p r + +newtype Arith p r a = Arith (ExceptT (Exception p r) (State (Context p r)) a) +instance Functor (Arith p r) +instance Applicative (Arith p r) +instance Monad (Arith p r) +evalArith :: Arith p r a -> Context p r -> Either (Exception p r) a + +data Signal + = Clamped + | DivisionByZero + | Inexact + | InvalidOperation + | Overflow + | Rounded + | Subnormal + | Underflow +raiseSignal :: Signal -> Decimal p r -> Arith p r (Decimal p r) + +type role Exception phantom phantom +data Exception p r +exceptionResult :: Exception p r -> Decimal p r 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 @@ -1,5 +1,5 @@ --- | The functions in this module implement conversions between 'Number' and +-- | The functions in this module implement conversions between 'Decimal' and -- 'String' as described in the /General Decimal Arithmetic Specification/. -- -- Because these functions are also used to implement 'Show' and 'Read' class @@ -28,11 +28,78 @@ import Numeric.Decimal.Number import Numeric.Decimal.Precision -import Numeric.Decimal.Rounding +{- $numeric-string-syntax + +(The following description is from the +/General Decimal Arithmetic Specification/.) + +Strings which are acceptable for conversion to the abstract representation of +numbers, or which might result from conversion from the abstract +representation to a string, are called /numeric strings/. + +A /numeric string/ is a character string that describes either a +/finite number/ or a /special value/. + +* If it describes a /finite number/, it includes one or more decimal + digits, with an optional decimal point. The decimal point may be embedded + in the digits, or may be prefixed or suffixed to them. The group of + digits (and optional point) thus constructed may have an optional sign + (“@+@” or “@-@”) which must come before any digits or decimal point. + + The string thus described may optionally be followed by an “@E@” + (indicating an exponential part), an optional sign, and an integer + following the sign that represents a power of ten that is to be + applied. The “@E@” may be in uppercase or lowercase. + +* If it describes a /special value/, it is one of the case-independent + names “@Infinity@”, “@Inf@”, “@NaN@”, or “@sNaN@” (where the first two + represent /infinity/ and the second two represent /quiet NaN/ and + /signaling NaN/ respectively). The name may be preceded by an optional + sign, as for finite numbers. If a NaN, the name may also be followed by + one or more digits, which encode any diagnostic information. + +No blanks or other white space characters are permitted in a numeric string. + +== Examples + +Some numeric strings are: + +> "0" -- zero +> "12" -- a whole number +> "-76" -- a signed whole number +> "12.70" -- some decimal places +> "+0.003" -- a plus sign is allowed, too +> "017." -- the same as 17 +> ".5" -- the same as 0.5 +> "4E+9" -- exponential notation +> "0.73e-7" -- exponential notation, negative power +> "Inf" -- the same as Infinity +> "-infinity" -- the same as -Inf +> "NaN" -- not-a-Number +> "NaN8275" -- diagnostic NaN + +== Notes + +1. A single period alone or with a sign is not a valid numeric string. +2. A sign alone is not a valid numeric string. +3. Significant (after the decimal point) and insignificant leading zeros + are permitted. +-} + +{- $setup +>>> :load Harness +>>> import Numeric.Decimal.Conversion +-} + +{- $doctest +prop> read' toNumber (show' toScientificString x) == x +prop> read' toNumber (show' toEngineeringString x) == x +-} + -- | Convert a number to a string, using scientific notation if an exponent is -- needed. -toScientificString :: Number p r -> ShowS +toScientificString :: Decimal p r -> ShowS toScientificString = showNumber exponential where exponential :: Exponent -> String -> Exponent -> ShowS @@ -40,9 +107,72 @@ showString ds . showExponent e exponential e ds _ = showString ds . showExponent e +{- $doctest-toScientificString + +>>> show' toScientificString $ fromRep $ N (0,123,0) +"123" + +>>> show' toScientificString $ fromRep $ N (1,123,0) +"-123" + +>>> show' toScientificString $ fromRep $ N (0,123,1) +"1.23E+3" + +>>> show' toScientificString $ fromRep $ N (0,123,3) +"1.23E+5" + +>>> show' toScientificString $ fromRep $ N (0,123,-1) +"12.3" + +>>> show' toScientificString $ fromRep $ N (0,123,-5) +"0.00123" + +>>> show' toScientificString $ fromRep $ N (0,123,-10) +"1.23E-8" + +>>> show' toScientificString $ fromRep $ N (1,123,-12) +"-1.23E-10" + +>>> show' toScientificString $ fromRep $ N (0,0,0) +"0" + +>>> show' toScientificString $ fromRep $ N (0,0,-2) +"0.00" + +>>> show' toScientificString $ fromRep $ N (0,0,2) +"0E+2" + +>>> show' toScientificString $ fromRep $ N (1,0,0) +"-0" + +>>> show' toScientificString $ fromRep $ N (0,5,-6) +"0.000005" + +>>> show' toScientificString $ fromRep $ N (0,50,-7) +"0.0000050" + +>>> show' toScientificString $ fromRep $ N (0,5,-7) +"5E-7" + +>>> show' toScientificString $ fromRep $ I 0 +"Infinity" + +>>> show' toScientificString $ fromRep $ I 1 +"-Infinity" + +>>> show' toScientificString $ fromRep $ Q (0,0) +"NaN" + +>>> show' toScientificString $ fromRep $ Q (0,123) +"NaN123" + +>>> show' toScientificString $ fromRep $ S (1,0) +"-sNaN" +-} + -- | Convert a number to a string, using engineering notation if an exponent -- is needed. -toEngineeringString :: Number p r -> ShowS +toEngineeringString :: Decimal p r -> ShowS toEngineeringString = showNumber exponential where exponential :: Exponent -> String -> Exponent -> ShowS @@ -70,8 +200,32 @@ shift _ e ds = showString ds . showExponent e +{- $doctest-toEngineeringString + +>>> show' toEngineeringString $ fromRep $ N (0,123,1) +"1.23E+3" + +>>> show' toEngineeringString $ fromRep $ N (0,123,3) +"123E+3" + +>>> show' toEngineeringString $ fromRep $ N (0,123,-10) +"12.3E-9" + +>>> show' toEngineeringString $ fromRep $ N (1,123,-12) +"-123E-12" + +>>> show' toEngineeringString $ fromRep $ N (0,7,-7) +"700E-9" + +>>> show' toEngineeringString $ fromRep $ N (0,7,1) +"70" + +>>> show' toEngineeringString $ fromRep $ N (0,0,1) +"0.00E+3" +-} + showNumber :: (Exponent -> String -> Exponent -> ShowS) - -> Number p r -> ShowS + -> Decimal p r -> ShowS showNumber exponential num = signStr . case num of Num { coefficient = c, exponent = e } | e <= 0 && ae >= -6 -> nonExponential @@ -114,21 +268,21 @@ -- | Convert a string to a number, as defined by its abstract representation. -- The string is expected to conform to the numeric string syntax described -- here. -toNumber :: (Precision p, Rounding r) => ReadP (Number p r) -toNumber = round <$> (parseSign flipSign <*> parseNumericString) +toNumber :: ReadP (Decimal PInfinite r) +toNumber = parseSign flipSign <*> parseNumericString where parseSign :: (a -> a) -> ReadP (a -> a) parseSign negate = char '-' *> pure negate <|> optional (char '+') *> pure id - parseNumericString :: ReadP (Number p r) + parseNumericString :: ReadP (Decimal p r) parseNumericString = parseNumericValue <|> parseNaN - parseNumericValue :: ReadP (Number p r) + parseNumericValue :: ReadP (Decimal p r) parseNumericValue = parseDecimalPart <*> option 0 parseExponentPart <|> parseInfinity - parseDecimalPart :: ReadP (Exponent -> Number p r) + parseDecimalPart :: ReadP (Exponent -> Decimal p r) parseDecimalPart = digitsWithPoint <|> digitsWithOptionalPoint where digitsWithPoint = do @@ -136,10 +290,9 @@ char '.' fracDigits <- many parseDigit return $ \e -> - Num { context = defaultContext - , sign = Pos + Num { sign = Pos , coefficient = readDigits (digits ++ fracDigits) - , exponent = e - fromIntegral (length fracDigits) + , exponent = e - fromIntegral (length fracDigits) } digitsWithOptionalPoint = fractionalDigits <|> wholeDigits @@ -148,18 +301,16 @@ char '.' fracDigits <- many1 parseDigit return $ \e -> - Num { context = defaultContext - , sign = Pos + Num { sign = Pos , coefficient = readDigits fracDigits - , exponent = e - fromIntegral (length fracDigits) + , exponent = e - fromIntegral (length fracDigits) } wholeDigits = do digits <- many1 parseDigit - return $ \e -> Num { context = defaultContext - , sign = Pos + return $ \e -> Num { sign = Pos , coefficient = readDigits digits - , exponent = e + , exponent = e } parseExponentPart :: ReadP Exponent @@ -167,25 +318,25 @@ parseString "E" parseSign negate <*> (readDigits <$> many1 parseDigit) - parseInfinity :: ReadP (Number p r) + parseInfinity :: ReadP (Decimal p r) parseInfinity = do parseString "Inf" optional $ parseString "inity" - return Inf { context = defaultContext, sign = Pos } + return Inf { sign = Pos } - parseNaN :: ReadP (Number p r) + parseNaN :: ReadP (Decimal p r) parseNaN = parseQNaN <|> parseSNaN - parseQNaN :: ReadP (Number p r) + parseQNaN :: ReadP (Decimal p r) parseQNaN = do p <- parseNaNPayload - return QNaN { context = defaultContext, sign = Pos, payload = p } + return QNaN { sign = Pos, payload = p } - parseSNaN :: ReadP (Number p r) + parseSNaN :: ReadP (Decimal p r) parseSNaN = do parseString "s" p <- parseNaNPayload - return SNaN { context = defaultContext, sign = Pos, payload = p } + return SNaN { sign = Pos, payload = p } parseNaNPayload :: ReadP Payload parseNaNPayload = do @@ -201,57 +352,74 @@ readDigits :: Num c => [Int] -> c readDigits = foldl' (\a b -> a * 10 + fromIntegral b) 0 --- $numeric-string-syntax --- --- (The following description is from the --- /General Decimal Arithmetic Specification/.) --- --- Strings which are acceptable for conversion to the abstract representation --- of numbers, or which might result from conversion from the abstract --- representation to a string, are called /numeric strings/. --- --- A /numeric string/ is a character string that describes either a /finite number/ or a /special value/. --- --- * If it describes a /finite number/, it includes one or more decimal --- digits, with an optional decimal point. The decimal point may be embedded --- in the digits, or may be prefixed or suffixed to them. The group of --- digits (and optional point) thus constructed may have an optional sign --- (“@+@” or “@-@”) which must come before any digits or decimal point. --- --- The string thus described may optionally be followed by an “@E@” --- (indicating an exponential part), an optional sign, and an integer --- following the sign that represents a power of ten that is to be --- applied. The “@E@” may be in uppercase or lowercase. --- --- * If it describes a /special value/, it is one of the case-independent --- names “@Infinity@”, “@Inf@”, “@NaN@”, or “@sNaN@” (where the first two --- represent /infinity/ and the second two represent /quiet NaN/ and --- /signaling NaN/ respectively). The name may be preceded by an optional --- sign, as for finite numbers. If a NaN, the name may also be followed by --- one or more digits, which encode any diagnostic information. --- --- No blanks or other white space characters are permitted in a numeric string. --- --- == Examples --- --- Some numeric strings are: --- --- > "0" -- zero --- > "12" -- a whole number --- > "-76" -- a signed whole number --- > "12.70" -- some decimal places --- > "+0.003" -- a plus sign is allowed, too --- > "017." -- the same as 17 --- > ".5" -- the same as 0.5 --- > "4E+9" -- exponential notation --- > "0.73e-7" -- exponential notation, negative power --- > "Inf" -- the same as Infinity --- > "-infinity" -- the same as -Inf --- > "NaN" -- not-a-Number --- > "NaN8275" -- diagnostic NaN --- --- == Notes --- --- 1. A single period alone or with a sign is not a valid numeric string. --- 2. A sign alone is not a valid numeric string. --- 3. Significant (after the decimal point) and insignificant leading zeros are permitted. +{- $doctest-toNumber + +>>> toRep $ read' toNumber "0" +N (0,0,0) + +>>> toRep $ read' toNumber "0.00" +N (0,0,-2) + +>>> toRep $ read' toNumber "123" +N (0,123,0) + +>>> toRep $ read' toNumber "-123" +N (1,123,0) + +>>> toRep $ read' toNumber "1.23E3" +N (0,123,1) + +>>> toRep $ read' toNumber "1.23E+3" +N (0,123,1) + +>>> toRep $ read' toNumber "12.3E+7" +N (0,123,6) + +>>> toRep $ read' toNumber "12.0" +N (0,120,-1) + +>>> toRep $ read' toNumber "12.3" +N (0,123,-1) + +>>> toRep $ read' toNumber "0.00123" +N (0,123,-5) + +>>> toRep $ read' toNumber "-1.23E-12" +N (1,123,-14) + +>>> toRep $ read' toNumber "1234.5E-4" +N (0,12345,-5) + +>>> toRep $ read' toNumber "-0" +N (1,0,0) + +>>> toRep $ read' toNumber "-0.00" +N (1,0,-2) + +>>> toRep $ read' toNumber "0E+7" +N (0,0,7) + +>>> toRep $ read' toNumber "-0E-7" +N (1,0,-7) + +>>> toRep $ read' toNumber "inf" +I 0 + +>>> toRep $ read' toNumber "+inFiniTy" +I 0 + +>>> toRep $ read' toNumber "-Infinity" +I 1 + +>>> toRep $ read' toNumber "NaN" +Q (0,0) + +>>> toRep $ read' toNumber "-NAN" +Q (1,0) + +>>> toRep $ read' toNumber "SNaN" +S (0,0) + +XXX toRep $ read' toNumber "Fred" +Q (0,0) +-} diff --git a/src/Numeric/Decimal/Conversion.hs-boot b/src/Numeric/Decimal/Conversion.hs-boot --- a/src/Numeric/Decimal/Conversion.hs-boot +++ b/src/Numeric/Decimal/Conversion.hs-boot @@ -8,8 +8,7 @@ import Text.ParserCombinators.ReadP (ReadP) import {-# SOURCE #-} Numeric.Decimal.Number -import Numeric.Decimal.Precision (Precision) -import {-# SOURCE #-} Numeric.Decimal.Rounding (Rounding) +import Numeric.Decimal.Precision -toScientificString :: Number p r -> ShowS -toNumber :: (Precision p, Rounding r) => ReadP (Number p r) +toScientificString :: Decimal p r -> ShowS +toNumber :: ReadP (Decimal PInfinite 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 @@ -10,7 +10,7 @@ , Exponent , Payload - , Number(..) + , Decimal(..) , zero , one , negativeOne @@ -20,7 +20,7 @@ , flipSign , cast - , excessDigits + , toBool , isPositive , isNegative @@ -28,36 +28,33 @@ , isZero , isNormal , isSubnormal - - , Context(..) - , TrapHandler - , defaultContext - , mergeContexts - - , Signal(..) - , raiseSignal ) where import Prelude hiding (exponent, round) -import Data.Bits (bit, complement, testBit, (.&.), (.|.)) +import Data.Char (isSpace) import Data.Coerce (coerce) -import Data.Monoid ((<>)) import Data.Ratio (numerator, denominator, (%)) import Numeric.Natural (Natural) import Text.ParserCombinators.ReadP (readP_to_S) import {-# SOURCE #-} Numeric.Decimal.Conversion +import {-# SOURCE #-} Numeric.Decimal.Arithmetic import Numeric.Decimal.Precision -import {-# SOURCE #-} Numeric.Decimal.Rounding +import Numeric.Decimal.Rounding import {-# SOURCE #-} qualified Numeric.Decimal.Operation as Op import qualified GHC.Real -data Sign = Pos | Neg - deriving (Eq, Enum, Show) +{- $setup +>>> :load Harness +-} +data Sign = Pos -- ^ Positive or non-negative + | Neg -- ^ Negative + deriving (Eq, Enum) + negateSign :: Sign -> Sign negateSign Pos = Neg negateSign Neg = Pos @@ -68,181 +65,365 @@ xorSigns Neg Pos = Neg xorSigns Neg Neg = Pos -signFactor :: Num a => Sign -> a -signFactor Pos = 1 -signFactor Neg = -1 - signFunc :: Num a => Sign -> a -> a signFunc Pos = id signFunc Neg = negate -type Coefficient = Natural -type Exponent = Integer +signMatch :: (Num a, Eq a) => a -> Sign +signMatch x = case signum x of + -1 -> Neg + _ -> Pos -type Payload = Coefficient +type Coefficient = Natural +type Exponent = Int +type Payload = Coefficient -- | A decimal floating point number with selectable precision and rounding -- algorithm -data Number p r - = Num { context :: Context p r - , sign :: Sign +data Decimal p r + = Num { sign :: Sign , coefficient :: Coefficient , exponent :: Exponent } - | Inf { context :: Context p r - , sign :: Sign + | Inf { sign :: Sign } - | QNaN { context :: Context p r - , sign :: Sign + | QNaN { sign :: Sign , payload :: Payload } - | SNaN { context :: Context p r - , sign :: Sign + | SNaN { sign :: Sign , payload :: Payload } -instance Precision p => Precision (Number p r) where - precision = precision . numberPrecision - where numberPrecision :: Number p r -> p - numberPrecision = undefined - -instance Show (Number p r) where +instance Show (Decimal p r) where showsPrec d n = showParen (d > 0 && isNegative n) $ toScientificString n -instance (Precision p, Rounding r) => Read (Number p r) where - readsPrec _ = readP_to_S toNumber +instance (Precision p, Rounding r) => Read (Decimal p r) where + readsPrec _ str = [ (cast n, s) + | (n, s) <- readParen False + (readP_to_S toNumber . dropWhile isSpace) str ] -instance (Precision p, Rounding r) => Eq (Number p r) where - x == y = case x `Op.compare` y of +{- $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)) +-} + +instance Precision p => Precision (Decimal p r) where + precision = precision . decimalPrecision + where decimalPrecision :: Decimal p r -> p + decimalPrecision = undefined + +evalOp :: (Precision p, Rounding r) => Arith p r (Decimal p r) -> Decimal p r +evalOp op = either exceptionResult id $ evalArith op newContext + +type GeneralDecimal = Decimal PInfinite RoundHalfEven + +instance Eq (Decimal p r) where + x == y = case evalOp (x `Op.compare` y) :: GeneralDecimal of Num { coefficient = 0 } -> True _ -> False -instance (Precision p, Rounding r) => Ord (Number p r) where - x `compare` y = case x `Op.compare` y of +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 - x < y = case x `Op.compare` y of + x < y = case evalOp (x `Op.compare` y) :: GeneralDecimal of Num { sign = Neg } -> True _ -> False - x <= y = case x `Op.compare` y of + x <= y = case evalOp (x `Op.compare` y) :: GeneralDecimal of Num { sign = Neg } -> True Num { coefficient = 0 } -> True _ -> False - x > y = case x `Op.compare` y of + x > y = case evalOp (x `Op.compare` y) :: GeneralDecimal of Num { coefficient = 0 } -> False Num { sign = Pos } -> True _ -> False - x >= y = case x `Op.compare` y of + x >= y = case evalOp (x `Op.compare` y) :: GeneralDecimal of Num { sign = Pos } -> True _ -> False - max nan@SNaN{} _ = nan - max _ nan@SNaN{} = nan - max nan@QNaN{} _ = nan - max _ nan@QNaN{} = nan - max x y - | x >= y = x - | otherwise = y + max x y = evalOp (Op.max x y) + min x y = evalOp (Op.min x y) - min nan@SNaN{} _ = nan - min _ nan@SNaN{} = nan - min nan@QNaN{} _ = nan - min _ nan@QNaN{} = nan - min x y - | x < y = x - | otherwise = 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) -instance (FinitePrecision p, Rounding r) => Enum (Number p r) where +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 +-} + +instance (Precision p, Rounding r) => Enum (Decimal p r) where + succ x = evalOp (x `Op.add` one) + pred x = evalOp (x `Op.subtract` one) + toEnum = fromIntegral - fromEnum = truncate -instance (Precision p, Rounding r) => Num (Number p r) where - (+) = Op.add - (-) = Op.subtract - (*) = Op.multiply - negate = Op.minus - abs = Op.abs + fromEnum Num { sign = s, coefficient = c, exponent = e } + | e >= 0 = signFunc s (fromIntegral c * 10^ e ) + | otherwise = signFunc s (fromIntegral $ c `quot` (10^(-e))) + fromEnum _ = 0 + enumFrom x = enumFromWith x one + enumFromThen x y = let i = y - x + in x : enumFromWith y i + enumFromTo x z = takeWhile (<= z) $ enumFromWith x one + enumFromThenTo x y z = let i = y - x + cmp | i < 0 = (>=) + | otherwise = (<=) + in takeWhile (`cmp` z) $ x : enumFromWith y i + +enumFromWith :: (Precision p, Rounding r) + => Decimal p r -> Decimal p r -> [Decimal p r] +enumFromWith x i = x : enumFromWith (x + i) i + +{- $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) + x * y = evalOp (x `Op.multiply` y) + + negate = evalOp . Op.minus + 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 - fromInteger x = Num { context = defaultContext - , sign = sx - , coefficient = fromInteger (abs x) - , exponent = 0 - } - where sx = case signum x of - -1 -> Neg - _ -> Pos + fromInteger x = cast + Num { sign = signMatch x + , coefficient = fromInteger (abs x) + , exponent = 0 + } -instance (Precision p, Rounding r) => Real (Number p r) where +{- $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 (signFactor s * fromIntegral c * 10^e) - | otherwise = (signFactor s * fromIntegral c) % 10^(-e) + | e >= 0 = fromInteger $ signFunc s (fromIntegral c * 10^e) + | otherwise = signFunc s (fromIntegral c) % 10^(-e) toRational n = signFunc (sign n) $ case n of Inf{} -> GHC.Real.infinity _ -> GHC.Real.notANumber -instance (FinitePrecision p, Rounding r) => Fractional (Number p r) where - (/) = Op.divide - fromRational r = fromInteger (numerator r) / fromInteger (denominator r) +instance (FinitePrecision p, Rounding r) => Fractional (Decimal p r) where + x / y = evalOp (x `Op.divide` y) + fromRational r = let n = fromInteger (numerator r) :: GeneralDecimal + d = fromInteger (denominator r) :: GeneralDecimal + in coerce n / coerce d -instance (FinitePrecision p, Rounding r) => RealFrac (Number p r) where +{- $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) - | otherwise = (signFactor s * fromIntegral c * 10^e, zero) - where n = signFactor s * fromIntegral q - f = x { coefficient = r, exponent = -(fromIntegral $ numDigits r) } + | otherwise = (signFunc s (fromIntegral c * 10^e), zero) + where n = signFunc s (fromIntegral q) + f = x { coefficient = r } (q, r) = c `quotRem` (10^(-e)) properFraction nan = (0, nan) --- | A 'Number' representing the value zero -zero :: Number p r -zero = Num { context = defaultContext - , sign = Pos +{- $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 an infinite series to maximum precision. +infiniteSeries :: (FinitePrecision p, Rounding r) + => [Decimal p r] -> Decimal p r +infiniteSeries = series zero + where series n (x:xs) + | n' == n = n' + | otherwise = series n' xs + where n' = n + 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) + x2 = x * x + +-- | Compute π to maximum precision using the arcsine series expansion. +seriesPi :: FinitePrecision p => Decimal p RoundHalfEven +seriesPi = 6 * arcsine oneHalf + +-- | 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 + +notyet :: String -> a +notyet = error . (++ ": not yet implemented") + +instance (FinitePrecision p, Rounding r) => Floating (Decimal p r) where + pi = castDown seriesPi + + exp = notyet "exp" + log = notyet "log" + + 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" + +{- $doctest-Floating +prop> realToFrac (pi :: ExtendedDecimal P16) == (pi :: Double) +-} + +instance (FinitePrecision p, Rounding r) => RealFloat (Decimal p r) where + floatRadix _ = 10 + floatDigits x = let Just p = precision x in p + floatRange _ = (minBound, maxBound) -- ? + + 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) + where special :: Sign -> Coefficient -> Integer + special s v = signFunc s (pp + fromIntegral v) + pp = 10 ^ floatDigits x :: Integer + + encodeFloat m n = x + where x | am >= pp = special + | otherwise = cast Num { sign = signMatch m + , coefficient = fromInteger am + , exponent = fromIntegral n + } + special + | n == maxBound = Inf { sign = signMatch m } + | n == minBound = QNaN { sign = signMatch m, payload = p } + | otherwise = SNaN { sign = signMatch m, payload = p } + where p = fromInteger (am - pp) + am = abs m :: Integer + pp = 10 ^ floatDigits x :: Integer + + isNaN x = case x of + QNaN{} -> True + SNaN{} -> True + _ -> False + + isInfinite x = case x of + Inf{} -> True + _ -> False + + isDenormalized = isSubnormal + + isNegativeZero x = case x of + Num { sign = Neg, coefficient = 0 } -> True + _ -> False + + isIEEE _ = True + +{- $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 +-} + +-- | A 'Decimal' representing the value zero +zero :: Decimal p r +zero = Num { sign = Pos , coefficient = 0 , exponent = 0 } --- | A 'Number' representing the value one -one :: Number p r +-- | A 'Decimal' representing the value ½ +oneHalf :: Decimal p r +oneHalf = zero { coefficient = 5, exponent = -1 } + +-- | A 'Decimal' representing the value one +one :: Decimal p r one = zero { coefficient = 1 } --- | A 'Number' representing the value negative one -negativeOne :: Number p r +-- | A 'Decimal' representing the value two +two :: Decimal p r +two = zero { coefficient = 2 } + +-- | A 'Decimal' representing the value negative one +negativeOne :: Decimal p r negativeOne = one { sign = Neg } --- | A 'Number' representing the value positive infinity -infinity :: Number p r -infinity = Inf { context = defaultContext, sign = Pos } +-- | A 'Decimal' representing the value positive infinity +infinity :: Decimal p r +infinity = Inf { sign = Pos } --- | A 'Number' representing undefined results -qNaN :: Number p r -qNaN = QNaN { context = defaultContext, sign = Pos, payload = 0 } +-- | A 'Decimal' representing undefined results +qNaN :: Decimal p r +qNaN = QNaN { sign = Pos, payload = 0 } --- | A signaling 'Number' representing undefined results -sNaN :: Number p r -sNaN = SNaN { context = defaultContext, sign = Pos, payload = 0 } +-- | A signaling 'Decimal' representing undefined results +sNaN :: Decimal p r +sNaN = SNaN { sign = Pos, payload = 0 } --- | Negate the given 'Number' by directly flipping its sign. -flipSign :: Number p r -> Number p r +-- | Negate the given 'Decimal' by directly flipping its sign. +flipSign :: Decimal p r -> Decimal p r flipSign n = n { sign = negateSign (sign n) } --- | Cast a 'Number' to another precision and/or rounding algorithm, +-- | Cast a 'Decimal' to another precision and/or rounding algorithm, -- immediately rounding if necessary to the new precision using the new -- algorithm. -cast :: (Precision p, Rounding r) => Number a b -> Number p r -cast = round . coerce +cast :: (Precision p, Rounding r) => Decimal a b -> Decimal p r +cast = evalOp . round . coerce +-- | Return the number of decimal digits of the argument. numDigits :: Coefficient -> Int numDigits x | x < 10 = 1 @@ -256,73 +437,71 @@ | x < 1000000000 = 9 | otherwise = 9 + numDigits (x `quot` 1000000000) -excessDigits :: Precision p => Number p r -> Maybe Int -excessDigits x@Num { coefficient = c } = precision x >>= excess - where excess p - | d > p = Just (d - p) - | otherwise = Nothing - where d = numDigits c -excessDigits _ = Nothing - maxCoefficient :: Precision p => p -> Maybe Coefficient maxCoefficient p = (\d -> 10 ^ d - 1) <$> precision p --- | Is the sign of the given 'Number' positive? -isPositive :: Number p r -> Bool +-- | Is the sign of the given 'Decimal' positive? +isPositive :: Decimal p r -> Bool isPositive n = case sign n of Pos -> True Neg -> False --- | Is the sign of the given 'Number' negative? -isNegative :: Number p r -> Bool +-- | Is the sign of the given 'Decimal' negative? +isNegative :: Decimal p r -> Bool isNegative n = case sign n of Neg -> True Pos -> False --- | Does the given 'Number' represent a finite value? -isFinite :: Number p r -> Bool +-- | Does the given 'Decimal' represent a finite value? +isFinite :: Decimal p r -> Bool isFinite Num{} = True isFinite _ = False --- | Does the given 'Number' represent the value zero? -isZero :: Number p r -> Bool +-- | Does the given 'Decimal' represent the value zero? +isZero :: Decimal p r -> Bool isZero Num { coefficient = 0 } = True isZero _ = False --- | Is the given 'Number' normal? -isNormal :: Precision p => Number p r -> Bool +-- | Is the given 'Decimal' normal? +isNormal :: Precision p => Decimal p r -> Bool isNormal n - | isFinite n && not (isZero n) && - maybe True (adjustedExponent n >=) (eMin n) = True - | otherwise = False + | isFinite n && not (isZero n) = maybe True (adjustedExponent n >=) (eMin n) + | otherwise = False --- | Is the given 'Number' subnormal? -isSubnormal :: Precision p => Number p r -> Bool +-- | Is the given 'Decimal' subnormal? +isSubnormal :: Precision p => Decimal p r -> Bool isSubnormal n - | isFinite n && not (isZero n) && - maybe False (adjustedExponent n <) (eMin n) = True - | otherwise = False + | isFinite n && not (isZero n) = maybe False (adjustedExponent n <) (eMin n) + | otherwise = False +-- | 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 + -- | Upper limit on the absolute value of the exponent -eLimit :: Precision p => Number p r -> Maybe Exponent -eLimit = eMax +eLimit :: Precision p => p -> Maybe Exponent +eLimit = eMax -- ? -- | Minimum value of the adjusted exponent -eMin :: Precision p => Number p r -> Maybe Exponent +eMin :: Precision p => p -> Maybe Exponent eMin n = (1 -) <$> eMax n -- | Maximum value of the adjusted exponent -eMax :: Precision p => Number p r -> Maybe 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 Natural + base = (10 *) . fromIntegral <$> mlength :: Maybe Coefficient -- | Minimum value of the exponent for subnormal results -eTiny :: Precision p => Number p r -> Maybe Exponent +eTiny :: Precision p => p -> Maybe Exponent eTiny n = (-) <$> eMin n <*> (fromIntegral . subtract 1 <$> precision n) -- | Range of permissible exponent values -eRange :: Precision p => Number p r -> Maybe (Exponent, Exponent) +eRange :: Precision p => Decimal p r -> Maybe (Exponent, Exponent) eRange n@Num { coefficient = c } = range <$> eLimit n where range :: Exponent -> (Exponent, Exponent) range lim = (-lim - clm1 + 1, lim - clm1) @@ -330,72 +509,8 @@ clm1 = fromIntegral (clength - 1) :: Exponent eRange _ = Nothing -adjustedExponent :: Number p r -> Exponent +adjustedExponent :: Decimal p r -> Exponent adjustedExponent Num { coefficient = c, exponent = e } = e + fromIntegral (clength - 1) where clength = numDigits c :: Int adjustedExponent _ = error "adjustedExponent: not a finite number" - -type TrapHandler p r = Signal -> Number p r -> Number p r - -data Context p r = Context { signalFlags :: Signals - , trapHandler :: TrapHandler p r - } - -instance Precision p => Precision (Context p r) where - precision = precision . contextPrecision - where contextPrecision :: Context p r -> p - contextPrecision = undefined - -defaultContext :: Context p r -defaultContext = Context mempty (const id) - -setSignal :: Signal -> Context p r -> Context p r -setSignal sig cxt = cxt { signalFlags = signalFlags cxt <> signal sig } - -modifyContext :: (Context p r -> Context p r) -> Number p r -> Number p r -modifyContext f n = n { context = f (context n) } - -mergeContexts :: Context p r -> Context p r -> Context p r -mergeContexts cxt1 cxt2 = - cxt1 { signalFlags = signalFlags cxt1 <> signalFlags cxt2 } - -data Signal - = Clamped - | DivisionByZero - | Inexact - | InvalidOperation - | Overflow - | Rounded - | Subnormal - | Underflow - deriving (Enum, Bounded, Show) - -newtype Signals = Signals Int - -instance Show Signals where - showsPrec d sigs = showParen (d > 10) $ - showString "signals " . showsPrec 11 (signalList sigs) - -instance Monoid Signals where - mempty = Signals 0 - Signals x `mappend` Signals y = Signals (x .|. y) - -signal :: Signal -> Signals -signal = Signals . bit . fromEnum - -unsignal :: Signal -> Signals -> Signals -unsignal sig (Signals ss) = Signals $ ss .&. complement (bit $ fromEnum sig) - -signals :: [Signal] -> Signals -signals = foldr (\s n -> signal s <> n) mempty - -signalList :: Signals -> [Signal] -signalList sigs = filter (testSignal sigs) [minBound..maxBound] - -testSignal :: Signals -> Signal -> Bool -testSignal (Signals ss) = testBit ss . fromEnum - -raiseSignal :: Signal -> Number p r -> Number p r -raiseSignal sig n = let n' = modifyContext (setSignal sig) n - in trapHandler (context n') sig n' 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 @@ -3,8 +3,38 @@ {-# LANGUAGE RoleAnnotations #-} module Numeric.Decimal.Number - ( Number + ( Sign(..) + , Decimal(..) + , Coefficient + , numDigits ) where -type role Number phantom phantom -data Number p r +import Numeric.Natural (Natural) + +import Numeric.Decimal.Precision + +data Sign = Pos | Neg +instance Eq Sign + +type Coefficient = Natural +type Exponent = Int +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) + +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 @@ -1,202 +1,1054 @@ --- | Eventually most or all of the arithmetic operations described in the --- /General Decimal Arithmetic Specification/ will be provided here. For now, --- the operations are mostly limited to those exposed through various class --- methods. --- --- It is not usually necessary to import this module. - -module Numeric.Decimal.Operation - ( abs - , add - , subtract - , multiply - , divide - , plus - , minus - , compare - ) where - -import Prelude hiding (abs, compare, exponent, round, subtract) -import qualified Prelude - -import Numeric.Decimal.Number -import Numeric.Decimal.Precision -import {-# SOURCE #-} Numeric.Decimal.Rounding - -invalidOperation :: Number p r -> Number p r -invalidOperation n = raiseSignal InvalidOperation qNaN { context = context n } - -toQNaN :: Number p r -> Number p r -toQNaN SNaN { context = t, sign = s, payload = p } = - QNaN { context = t, sign = s, payload = p } -toQNaN n@QNaN{} = n -toQNaN n = qNaN { context = context n, sign = sign n } - -toQNaN2 :: Number p r -> Number p r -> Number p r -toQNaN2 nan@SNaN{} _ = toQNaN nan -toQNaN2 _ nan@SNaN{} = toQNaN nan -toQNaN2 nan@QNaN{} _ = nan -toQNaN2 _ nan@QNaN{} = nan -toQNaN2 n _ = toQNaN n - --- | Add two operands. -add :: (Precision p, Rounding r) => Number p r -> Number p r -> Number p r -add Num { context = xt, sign = xs, coefficient = xc, exponent = xe } - Num { context = yt, sign = ys, coefficient = yc, exponent = ye } = round rn - - where rn = Num { context = rt, sign = rs, coefficient = rc, exponent = re } - rt = mergeContexts xt yt - rs | rc /= 0 = if xac > yac then xs else ys - | xs == Neg && ys == Neg = Neg - | xs /= ys && isRoundFloor rn = Neg - | otherwise = Pos - rc | xs == ys = xac + yac - | xac > yac = xac - yac - | otherwise = yac - xac - re = 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 { context = xt, sign = xs } Inf { context = yt, sign = ys } - | xs == ys = inf { context = mergeContexts xt yt } - | otherwise = invalidOperation inf { context = mergeContexts xt yt } -add inf@Inf{} Num{} = inf -add Num{} inf@Inf{} = inf -add x y = toQNaN2 x y - --- | Subtract the second operand from the first. -subtract :: (Precision p, Rounding r) => Number p r -> Number p r -> Number p r -subtract x = add x . flipSign - --- | Unary minus (negation) -minus :: (Precision p, Rounding r) => Number p r -> Number p r -minus x = zero { exponent = exponent x } `subtract` x - --- | Unary plus -plus :: (Precision p, Rounding r) => Number p r -> Number p r -plus x = zero { exponent = exponent x } `add` x - --- | Multiply two operands. -multiply :: (Precision p, Rounding r) => Number p r -> Number p r -> Number p r -multiply Num { context = xt, sign = xs, coefficient = xc, exponent = xe } - Num { context = yt, sign = ys, coefficient = yc, exponent = ye } = - round rn - - where rn = Num { context = rt, sign = rs, coefficient = rc, exponent = re } - rt = mergeContexts xt yt - rs = xorSigns xs ys - rc = xc * yc - re = xe + ye - -multiply Inf { context = xt, sign = xs } Inf { context = yt, sign = ys } = - Inf { context = mergeContexts xt yt, sign = xorSigns xs ys } -multiply Inf { context = xt, sign = xs } Num { context = yt, sign = ys } = - Inf { context = mergeContexts xt yt, sign = xorSigns xs ys } -multiply Num { context = xt, sign = xs } Inf { context = yt, sign = ys } = - Inf { context = mergeContexts xt yt, sign = xorSigns xs ys } -multiply x y = toQNaN2 x y - --- | Divide the first dividend operand by the second divisor using long division. -divide :: (FinitePrecision p, Rounding r) - => Number p r -> Number p r -> Number p r -divide dividend@Num{ sign = xs } Num { coefficient = 0, sign = ys } - | isZero dividend = invalidOperation qNaN - | otherwise = raiseSignal DivisionByZero - infinity { sign = xorSigns xs ys } -divide Num { context = xt, sign = xs, coefficient = xc, exponent = xe } - Num { context = yt, sign = ys, coefficient = yc, exponent = ye } = - result - - where rn = Num { context = rt, sign = rs, coefficient = rc, exponent = re } - rt = mergeContexts xt yt - rs = xorSigns xs ys - (rc, rem, dv, adjust) = longDivision xc yc p - re = xe - (ye + adjust) - Just p = precision rn - result - | rem == 0 = rn - | otherwise = round $ 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 { context = xt, sign = xs } Num { context = yt, sign = ys } = - Inf { context = mergeContexts xt yt, sign = xorSigns xs ys } -divide Num { context = xt, sign = xs } Inf { context = yt, sign = ys } = - zero { context = mergeContexts xt yt, sign = xorSigns xs ys } -divide x y = toQNaN2 x y - -type Dividend = Coefficient -type Divisor = Coefficient -type Quotient = Coefficient -type Remainder = Coefficient - -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 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 = step3 0 - - 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 = (,,,) - --- | 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. -abs :: (Precision p, Rounding r) => Number p r -> Number p r -abs x - | isNegative x = minus x - | otherwise = plus x - --- | Compare the values of two operands numerically, 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) => Number p r -> Number p r -> Number p r -compare x@Num{} y@Num{} = (nzp $ xn `subtract` yn) { context = rt } - - where (xn, yn) | sign x /= sign y = (nzp x, nzp y) - | otherwise = (x, y) - - rt = mergeContexts (context x) (context y) - - nzp :: Number p r -> Number p r - nzp Num { context = t, sign = s, coefficient = c } - | c == 0 = zero { context = t } - | s == Pos = one { context = t } - | otherwise = negativeOne { context = t } - nzp Inf { context = t, sign = s } - | s == Pos = one { context = t } - | otherwise = negativeOne { context = t } - nzp n = toQNaN n - -compare Inf { context = xt, sign = xs } Inf { context = yt, sign = ys } - | xs == ys = zero { context = rt } - | xs == Neg = negativeOne { context = rt } - | otherwise = one { context = rt } - where rt = mergeContexts xt yt -compare Inf { context = xt, sign = xs } Num { context = yt } - | xs == Neg = negativeOne { context = rt } - | otherwise = one { context = rt } - where rt = mergeContexts xt yt -compare Num { context = xt } Inf { context = yt, sign = ys } - | ys == Pos = negativeOne { context = rt } - | otherwise = one { context = rt } - where rt = mergeContexts xt yt -compare nan@SNaN{} _ = invalidOperation nan -compare _ nan@SNaN{} = invalidOperation nan -compare x y = toQNaN2 x y +{- | Eventually most or all of the arithmetic operations described in the +/General Decimal Arithmetic Specification/ will be provided here. For now, the +operations are mostly limited to those exposed through various class methods. + +It is suggested to import this module qualified to avoid "Prelude" name +clashes: + +> import qualified Numeric.Decimal.Operation as Op + +Note that it is not usually necessary to import this module unless you want to +use operations unavailable through class methods, or you need precise control +over the handling of exceptional conditions. +-} +module Numeric.Decimal.Operation + ( -- * 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(..), Subclass(..) + -- 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, compare, exponent, isInfinite, isNaN, max, min, + round, subtract) +import qualified Prelude + +import Data.Coerce (coerce) + +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 = round -- ... +-- | 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 + +-- $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 } = + return Inf { sign = xorSigns xs ys } +multiply Num { sign = xs } Inf { sign = ys } = + return Inf { sign = xorSigns xs ys } +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 +-} + +-- | '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 = round $ 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 = (,,,) + +-- | '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 +-} + +-- | '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 } + +-- | '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 +-} + +-- $miscellaneous-operations +-- +-- This section describes miscellaneous operations on decimal numbers, +-- including non-numeric comparisons, sign and other manipulations, and +-- logical operations. + +-- | '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 -> Class (sign n) ZeroClass + | Number.isSubnormal n -> Class (sign n) SubnormalClass + | otherwise -> Class (sign n) NormalClass + Inf {} -> Class (sign n) InfinityClass + QNaN{} -> Class Pos NaNClass + SNaN{} -> Class Neg NaNClass + +data Class = Class Sign Subclass deriving Eq + +data Subclass = ZeroClass -- ^ Zero + | NormalClass -- ^ Normal finite number + | SubnormalClass -- ^ Subnormal finite number + | InfinityClass -- ^ Infinity + | NaNClass -- ^ Not a number (quiet or signaling) + deriving Eq + +instance Show Class where + show c = case c of + Class Pos s@NaNClass -> showSubclass s + Class Neg s@NaNClass -> 's' : showSubclass s + Class Pos s -> '+' : showSubclass s + Class Neg s -> '-' : showSubclass s + + where showSubclass s = case s of + ZeroClass -> "Zero" + NormalClass -> "Normal" + SubnormalClass -> "Subnormal" + InfinityClass -> "Infinity" + NaNClass -> "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 (Decimal p r) +isCanonical _ = return one + +{- $doctest-isCanonical +>>> 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 (Decimal p r) +isFinite n = return $ case n of + Num{} -> one + _ -> zero + +{- $doctest-isFinite +>>> op1 Op.isFinite "2.50" +1 + +>>> op1 Op.isFinite "-0.3" +1 + +>>> op1 Op.isFinite "0" +1 + +>>> op1 Op.isFinite "Inf" +0 + +>>> 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 (Decimal p r) +isInfinite n = return $ case n of + Inf{} -> one + _ -> zero + +{- $doctest-isInfinite +>>> op1 Op.isInfinite "2.50" +0 + +>>> op1 Op.isInfinite "-Inf" +1 + +>>> 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 (Decimal p r) +isNaN n = return $ case n of + QNaN{} -> one + SNaN{} -> one + _ -> zero + +{- $doctest-isNaN +>>> op1 Op.isNaN "2.50" +0 + +>>> op1 Op.isNaN "NaN" +1 + +>>> 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 (Decimal p r) +isNormal n = return $ case n of + _ | Number.isNormal n -> one + | otherwise -> zero + +{- $doctest-isNormal +>>> op1 Op.isNormal "2.50" +1 + +>>> op1 Op.isNormal "0.1E-999" +0 + +>>> op1 Op.isNormal "0.00" +0 + +>>> op1 Op.isNormal "-Inf" +0 + +>>> 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 (Decimal p r) +isQNaN n = return $ case n of + QNaN{} -> one + _ -> zero + +{- $doctest-isQNaN +>>> op1 Op.isQNaN "2.50" +0 + +>>> op1 Op.isQNaN "NaN" +1 + +>>> 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 (Decimal p r) +isSigned n = return $ case sign n of + Neg -> one + Pos -> zero + +{- $doctest-isSigned +>>> op1 Op.isSigned "2.50" +0 + +>>> op1 Op.isSigned "-12" +1 + +>>> 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 (Decimal p r) +isSNaN n = return $ case n of + SNaN{} -> one + _ -> zero + +{- $doctest-isSNaN +>>> op1 Op.isSNaN "2.50" +0 + +>>> op1 Op.isSNaN "NaN" +0 + +>>> 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 (Decimal p r) +isSubnormal n = return $ case n of + _ | Number.isSubnormal n -> one + | otherwise -> zero + +{- $doctest-isSubnormal +>>> op1 Op.isSubnormal "2.50" +0 + +>>> op1 Op.isSubnormal "0.1E-999" +1 + +>>> op1 Op.isSubnormal "0.00" +0 + +>>> op1 Op.isSubnormal "-Inf" +0 + +>>> 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 (Decimal p r) +isZero n = return $ case n of + _ | Number.isZero n -> one + | otherwise -> zero + +{- $doctest-isZero +>>> op1 Op.isZero "0" +1 + +>>> op1 Op.isZero "2.50" +0 + +>>> op1 Op.isZero "-0E+2" +1 +-} + +-- | '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 (Decimal p r) +sameQuantum Num { exponent = e1 } Num { exponent = e2 } + | e1 == e2 = return one + | otherwise = return zero +sameQuantum Inf {} Inf {} = return one +sameQuantum QNaN{} QNaN{} = return one +sameQuantum SNaN{} SNaN{} = return one +sameQuantum QNaN{} SNaN{} = return one +sameQuantum SNaN{} QNaN{} = return one +sameQuantum _ _ = return zero + +{- $doctest-sameQuantum +>>> op2 Op.sameQuantum "2.17" "0.001" +0 + +>>> op2 Op.sameQuantum "2.17" "0.01" +1 + +>>> op2 Op.sameQuantum "2.17" "0.1" +0 + +>>> op2 Op.sameQuantum "2.17" "1" +0 + +>>> op2 Op.sameQuantum "Inf" "-Inf" +1 + +>>> op2 Op.sameQuantum "NaN" "NaN" +1 +-} 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 @@ -8,19 +8,32 @@ , minus , abs , compare + , min + , max ) where -import Prelude hiding (abs, compare, subtract) +import Prelude hiding (abs, compare, max, min, subtract) +import {-# SOURCE #-} Numeric.Decimal.Arithmetic import {-# SOURCE #-} Numeric.Decimal.Number import Numeric.Decimal.Precision -import {-# SOURCE #-} Numeric.Decimal.Rounding +import Numeric.Decimal.Rounding -add :: (Precision p, Rounding r) => Number p r -> Number p r -> Number p r -subtract :: (Precision p, Rounding r) => Number p r -> Number p r -> Number p r -multiply :: (Precision p, Rounding r) => Number p r -> Number p r -> Number p r -divide :: (FinitePrecision p, Rounding r) => - Number p r -> Number p r -> Number p r -minus :: (Precision p, Rounding r) => Number p r -> Number p r -abs :: (Precision p, Rounding r) => Number p r -> Number p r -compare :: (Precision p, Rounding r) => Number p r -> Number p r -> Number p r +add :: (Precision p, Rounding r) + => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +subtract :: (Precision p, Rounding r) + => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +multiply :: (Precision p, Rounding r) + => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +divide :: (FinitePrecision p, Rounding r) + => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) +minus :: (Precision p, Rounding r) + => 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) 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 @@ -11,18 +11,19 @@ , P75, P100, P150, P200, P250, P300, P400, P500, P1000, P2000 - , PPlus1, PTimes2 + , PPlus1, PPlus2, PPlus3, PPlus4, PPlus5, PPlus6, PPlus7, PPlus8, PPlus9 + , PTimes2, PTimes10 , PInfinite ) where --- | Precision indicates the maximum number of significant digits a number may --- have. +-- | 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. precision :: p -> Maybe Int --- | A subclass of precisions which are finite +-- | A subclass of precisions that are finite class Precision p => FinitePrecision p -- | A precision of unlimited significant digits @@ -36,25 +37,43 @@ precision _ = Just 1 instance FinitePrecision P1 --- | A precision of (@p@ + 1) significant digits -data PPlus1 p -instance Precision p => Precision (PPlus1 p) where - precision pp = (+ 1) <$> precision (minus1 pp) - where minus1 :: PPlus1 p -> p - minus1 = undefined -instance FinitePrecision p => FinitePrecision (PPlus1 p) +-- | A precision of (@p@ + @q@) significant digits +data PPlus p q +instance (Precision p, Precision q) => Precision (PPlus p q) where + precision pq = (+) <$> precision (p pq) <*> precision (q pq) + where p :: PPlus p q -> p + p = undefined + q :: PPlus p q -> q + q = undefined +instance (FinitePrecision p, FinitePrecision q) => FinitePrecision (PPlus p q) --- | A precision of (@p@ × 2) significant digits -data PTimes2 p -instance Precision p => Precision (PTimes2 p) where - precision pp = (* 2) <$> precision (div2 pp) - where div2 :: PTimes2 p -> p - div2 = undefined -instance FinitePrecision p => FinitePrecision (PTimes2 p) +-- | A precision of (@p@ × @q@) significant digits +data PTimes p q +instance (Precision p, Precision q) => Precision (PTimes p q) where + precision pq = (*) <$> precision (p pq) <*> precision (q pq) + where p :: PTimes p q -> p + p = undefined + q :: PTimes p q -> q + q = undefined +instance (FinitePrecision p, FinitePrecision q) => FinitePrecision (PTimes p q) +type PPlus1 p = PPlus p P1 -- ^ A precision of (@p@ + 1) significant digits +type PPlus2 p = PPlus p P2 -- ^ A precision of (@p@ + 2) significant digits +type PPlus3 p = PPlus p P3 -- ^ A precision of (@p@ + 3) significant digits +type PPlus4 p = PPlus p P4 -- ^ A precision of (@p@ + 4) significant digits +type PPlus5 p = PPlus p P5 -- ^ A precision of (@p@ + 5) significant digits +type PPlus6 p = PPlus p P6 -- ^ A precision of (@p@ + 6) significant digits +type PPlus7 p = PPlus p P7 -- ^ A precision of (@p@ + 7) significant digits +type PPlus8 p = PPlus p P8 -- ^ A precision of (@p@ + 8) significant digits +type PPlus9 p = PPlus p P9 -- ^ A precision of (@p@ + 9) significant digits + +type PTimes2 = PTimes P2 +type PTimes10 = PTimes P10 + -- | A precision of 2 significant digits -type P2 = PTimes2 P1 ; type P3 = PPlus1 P2 --- ^ A precision of 3 significant digits +type P2 = PPlus1 P1 +-- | A precision of 3 significant digits +type P3 = PPlus1 P2 -- | Et cetera type P4 = PTimes2 P2 ; type P5 = PPlus1 P4 @@ -82,15 +101,13 @@ type P48 = PTimes2 P24; type P49 = PPlus1 P48 type P50 = PTimes2 P25 -type P62 = PTimes2 P31 -type P74 = PTimes2 P37; type P75 = PPlus1 P74 +type P75 = PPlus5 (PTimes10 P7) type P100 = PTimes2 P50 -type P124 = PTimes2 P62; type P125 = PPlus1 P124 type P150 = PTimes2 P75 type P200 = PTimes2 P100 -type P250 = PTimes2 P125 +type P250 = PTimes10 P25 type P300 = PTimes2 P150 type P400 = PTimes2 P200 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,6 +1,7 @@ module Numeric.Decimal.Rounding - ( Rounding(..) + ( RoundingAlgorithm(..) + , Rounding(..) , RoundDown , RoundHalfUp @@ -15,22 +16,32 @@ import Prelude hiding (exponent) -import Numeric.Decimal.Number -import Numeric.Decimal.Precision +import {-# SOURCE #-} Numeric.Decimal.Number +import Numeric.Decimal.Precision +import {-# SOURCE #-} Numeric.Decimal.Arithmetic +data RoundingAlgorithm = RoundDown + | RoundHalfUp + | RoundHalfEven + | RoundCeiling + | RoundFloor + | RoundHalfDown + | RoundUp + | Round05Up + deriving (Eq, Enum) + -- | A rounding algorithm to use when the result of an arithmetic operation -- exceeds the precision of the result type class Rounding r where - round :: Precision p => Number p r -> Number p r - - isRoundFloor :: Number p r -> Bool - isRoundFloor _ = False + rounding :: r -> RoundingAlgorithm + round :: Precision p => Decimal p r -> Arith p r (Decimal p r) -- Required... -- | Round toward 0 (truncate) data RoundDown instance Rounding RoundDown where + rounding _ = RoundDown round = roundDown -- | If the discarded digits represent greater than or equal to half (0.5) of @@ -38,6 +49,7 @@ -- up. If they represent less than half, the value is rounded down. data RoundHalfUp instance Rounding RoundHalfUp where + rounding _ = RoundHalfUp round = roundHalfUp -- | If the discarded digits represent greater than half (0.5) of the value of @@ -46,18 +58,20 @@ -- exactly half, the value is rounded to make its rightmost digit even. data RoundHalfEven instance Rounding RoundHalfEven where + rounding _ = RoundHalfEven round = roundHalfEven -- | Round toward +∞ data RoundCeiling instance Rounding RoundCeiling where + rounding _ = RoundCeiling round = roundCeiling -- | Round toward −∞ data RoundFloor instance Rounding RoundFloor where + rounding _ = RoundFloor round = roundFloor - isRoundFloor _ = True -- Optional... @@ -66,107 +80,115 @@ -- represent less than half or exactly half, the value is rounded down. data RoundHalfDown instance Rounding RoundHalfDown where + rounding _ = RoundHalfDown round = roundHalfDown -- | Round away from 0 data RoundUp instance Rounding RoundUp where + rounding _ = RoundUp round = roundUp -- | Round zero or five away from 0 data Round05Up instance Rounding Round05Up where + rounding _ = Round05Up round = round05Up -- Implementations +excessDigits :: Precision p => Decimal p r -> Arith p r (Maybe Int) +excessDigits n@Num { coefficient = c } = result + where result = return (precision n >>= excess) + d = numDigits c + excess p + | d > p = Just (d - p) + | otherwise = Nothing +excessDigits _ = return Nothing + rounded :: (Coefficient -> Coefficient -> Coefficient -> - Number p r -> Number p r -> Number p r) - -> Int -> Number p r -> Number p r -rounded f d n = raiseSignal Rounded rounded' + Decimal p r -> Decimal p r -> Decimal p r) + -> Int -> Decimal p r -> Arith p r (Decimal p r) +rounded f d n = raiseSignal Rounded =<< rounded' n' where rounded' - | r /= 0 = raiseSignal Inexact n' - | otherwise = n' + | r /= 0 = raiseSignal Inexact + | otherwise = return p = 10 ^ d (q, r) = coefficient n `quotRem` p n' = f (p `quot` 2) q r down up - down = n { coefficient = q - , exponent = exponent n + fromIntegral d - } - up = n { coefficient = q + 1 - , exponent = exponent n + fromIntegral d - } + down = n { coefficient = q , exponent = exponent n + fromIntegral d } + up = n { coefficient = q + 1, exponent = exponent n + fromIntegral d } -roundDown :: Precision p => Number p r -> Number p r -roundDown n = roundDown' (excessDigits n) - where roundDown' Nothing = n +roundDown :: Precision p => Decimal p r -> Arith p r (Decimal p r) +roundDown n = excessDigits n >>= roundDown' + where roundDown' Nothing = return n roundDown' (Just d) = rounded choice d n choice _h _q _r down _up = down -roundHalfUp :: Precision p => Number p r -> Number p r -roundHalfUp n = roundHalfUp' (excessDigits n) - where roundHalfUp' Nothing = n +roundHalfUp :: Precision p => Decimal p r -> Arith p r (Decimal p r) +roundHalfUp n = excessDigits n >>= roundHalfUp' + where roundHalfUp' Nothing = return n roundHalfUp' (Just d) = rounded choice d n choice h _q r down up - | r >= h = roundHalfUp up + | r >= h = up | otherwise = down -roundHalfEven :: Precision p => Number p r -> Number p r -roundHalfEven n = roundHalfEven' (excessDigits n) - where roundHalfEven' Nothing = n +roundHalfEven :: Precision p => Decimal p r -> Arith p r (Decimal p r) +roundHalfEven n = excessDigits n >>= roundHalfEven' + where roundHalfEven' Nothing = return n roundHalfEven' (Just d) = rounded choice d n choice h q r down up = case r `Prelude.compare` h of LT -> down - GT -> roundHalfEven up + GT -> up EQ | even q -> down - | otherwise -> roundHalfEven up + | otherwise -> up -roundCeiling :: Precision p => Number p r -> Number p r -roundCeiling n = roundCeiling' (excessDigits n) - where roundCeiling' Nothing = n +roundCeiling :: Precision p => Decimal p r -> Arith p r (Decimal p r) +roundCeiling n = excessDigits n >>= roundCeiling' + where roundCeiling' Nothing = return n roundCeiling' (Just d) = rounded choice d n choice _h _q r down up | r == 0 || sign n == Neg = down - | otherwise = roundCeiling up + | otherwise = up -roundFloor :: Precision p => Number p r -> Number p r -roundFloor n = roundFloor' (excessDigits n) - where roundFloor' Nothing = n +roundFloor :: Precision p => Decimal p r -> Arith p r (Decimal p r) +roundFloor n = excessDigits n >>= roundFloor' + where roundFloor' Nothing = return n roundFloor' (Just d) = rounded choice d n choice _h _q r down up | r == 0 || sign n == Pos = down - | otherwise = roundFloor up + | otherwise = up -roundHalfDown :: Precision p => Number p r -> Number p r -roundHalfDown n = roundHalfDown' (excessDigits n) - where roundHalfDown' Nothing = n +roundHalfDown :: Precision p => Decimal p r -> Arith p r (Decimal p r) +roundHalfDown n = excessDigits n >>= roundHalfDown' + where roundHalfDown' Nothing = return n roundHalfDown' (Just d) = rounded choice d n choice h _q r down up - | r > h = roundHalfDown up + | r > h = up | otherwise = down -roundUp :: Precision p => Number p r -> Number p r -roundUp n = roundUp' (excessDigits n) - where roundUp' Nothing = n +roundUp :: Precision p => Decimal p r -> Arith p r (Decimal p r) +roundUp n = excessDigits n >>= roundUp' + where roundUp' Nothing = return n roundUp' (Just d) = rounded choice d n choice _h _q r down up | r == 0 = down - | otherwise = roundUp up + | otherwise = up -round05Up :: Precision p => Number p r -> Number p r -round05Up n = round05Up' (excessDigits n) - where round05Up' Nothing = n +round05Up :: Precision p => Decimal p r -> Arith p r (Decimal p r) +round05Up n = excessDigits n >>= round05Up' + where round05Up' Nothing = return n round05Up' (Just d) = rounded choice d n choice _h q r down up | r == 0 = down - | d == 0 || d == 5 = round05Up up -- overflow -> roundDown? + | d == 0 || d == 5 = up -- XXX overflow -> roundDown? | otherwise = down where d = q `rem` 10 diff --git a/src/Numeric/Decimal/Rounding.hs-boot b/src/Numeric/Decimal/Rounding.hs-boot deleted file mode 100644 --- a/src/Numeric/Decimal/Rounding.hs-boot +++ /dev/null @@ -1,14 +0,0 @@ --- -*- Haskell -*- - -module Numeric.Decimal.Rounding - ( Rounding(..) - ) where - -import {-# SOURCE #-} Numeric.Decimal.Number (Number) -import Numeric.Decimal.Precision (Precision) - -class Rounding r where - round :: Precision p => Number p r -> Number p r - - isRoundFloor :: Number p r -> Bool - isRoundFloor _ = False diff --git a/test/Spec.hs b/test/Spec.hs deleted file mode 100644 --- a/test/Spec.hs +++ /dev/null @@ -1,30 +0,0 @@ - -import Numeric.Decimal -import Test.QuickCheck - -main :: IO () -main = putStrLn "Test suite not yet implemented" - -infinity :: (Precision p, Rounding r) => Number p r -infinity = read "Infinity" - -instance (Precision p, Rounding r) => Arbitrary (Number p r) where - arbitrary = frequency [(85, genNum), (10, genInf)] - -genNum :: (Precision p, Rounding r) => Gen (Number 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 (Number p r) -genInf = do - s <- elements [-1, 1] - return (s * infinity) - -genNaN :: (Precision p, Rounding r) => Gen (Number 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/doctests.hs b/test/doctests.hs new file mode 100644 --- /dev/null +++ b/test/doctests.hs @@ -0,0 +1,4 @@ + +import Test.DocTest + +main = doctest ["-isrc:test", "src/Numeric/Decimal.hs"]