decimal-arithmetic 0.4.0.0 → 0.5.0.0

raw patch · 21 files changed

+3316/−2254 lines, 21 filesdep +binarydep +binary-bitsdep +deepseqdep ~basedep ~mtlPVP ok

version bump matches the API change (PVP)

Dependencies added: binary, binary-bits, deepseq, hspec

Dependency ranges changed: base, mtl

API changes (from Hackage documentation)

Files

• LICENSE +1/−1
• TODO +0/−1
• decimal-arithmetic.cabal +27/−6
• src/Numeric/Decimal.hs +66/−29
• src/Numeric/Decimal/Arithmetic.hs +8/−9
• src/Numeric/Decimal/Conversion.hs +16/−21
• src/Numeric/Decimal/Encoding.hs +330/−0
• src/Numeric/Decimal/Exception.hs +242/−0
• src/Numeric/Decimal/Exception.hs-boot +11/−0
• src/Numeric/Decimal/Number.hs +329/−190
• src/Numeric/Decimal/Number.hs-boot +13/−19
• src/Numeric/Decimal/Operation.hs +1380/−1947
• src/Numeric/Decimal/Operation.hs-boot +6/−6
• src/Numeric/Decimal/Precision.hs +17/−1
• src/Numeric/Decimal/Rounding.hs +18/−12
• stack.yaml +12/−12
• test/Arbitrary.hs +31/−0
• test/Numeric/Decimal/EncodingSpec.hs +25/−0
• test/Numeric/Decimal/NumberSpec.hs +165/−0
• test/Numeric/Decimal/OperationSpec.hs +618/−0
• test/Spec.hs +1/−0

@@ -1,4 +1,4 @@-Copyright (c) 2016, Robert Leslie+Copyright (c) 2016-2017, Robert Leslie  All rights reserved. 

@@ -1,5 +1,4 @@ -*- Outline -*-  * To Do-** instance Floating (Decimal p r) ** instance PrintfArg (Decimal p r)

@@ -1,6 +1,6 @@  name:                decimal-arithmetic-version:             0.4.0.0+version:             0.5.0.0  synopsis:            An implementation of the General Decimal Arithmetic                      Specification@@ -16,7 +16,7 @@ license:             BSD3 license-file:        LICENSE -copyright:           © 2016 Robert Leslie+copyright:           © 2016–2017 Robert Leslie author:              Rob Leslie <rob@mars.org> maintainer:          Rob Leslie <rob@mars.org> @@ -38,20 +38,41 @@   hs-source-dirs:      src    exposed-modules:     Numeric.Decimal-                       Numeric.Decimal.Conversion                        Numeric.Decimal.Arithmetic+                       Numeric.Decimal.Conversion+                       Numeric.Decimal.Encoding                        Numeric.Decimal.Operation-  other-modules:       Numeric.Decimal.Number+  other-modules:       Numeric.Decimal.Exception+                       Numeric.Decimal.Number                        Numeric.Decimal.Precision                        Numeric.Decimal.Rounding -  build-depends:       base >= 4.7 && < 5-                     , mtl+  build-depends:       base >= 4.8 && < 5+                     , binary >= 0.8 && < 0.9+                     , binary-bits >= 0.5 && < 0.6+                     , deepseq >= 1.4 && < 1.5+                     , mtl >= 2.2 && < 2.3   default-language:    Haskell2010   default-extensions:  Trustworthy   other-extensions:    FlexibleInstances                        MultiParamTypeClasses                        RoleAnnotations++test-suite spec+  type:                exitcode-stdio-1.0+  hs-source-dirs:      test+  main-is:             Spec.hs+  other-modules:       Arbitrary+                       Numeric.Decimal.EncodingSpec+                       Numeric.Decimal.NumberSpec+                       Numeric.Decimal.OperationSpec+  build-depends:       base+                     , binary+                     , decimal-arithmetic+                     , hspec+                     , QuickCheck+  ghc-options:         -threaded -rtsopts -with-rtsopts=-N+  default-language:    Haskell2010  test-suite doctests   type:                exitcode-stdio-1.0

@@ -1,7 +1,7 @@ {-| Module      : Numeric.Decimal Description : General arbitrary-precision decimal floating-point number type-Copyright   : © 2016 Robert Leslie+Copyright   : © 2016–2017 Robert Leslie License     : BSD3 Maintainer  : rob@mars.org Stability   : experimental@@ -38,20 +38,49 @@ representation via 'Show' and 'Read' instances. Note that there may be multiple representations of values that are numerically equal (e.g. 1 and 1.00) which are preserved by this conversion.++Some decimal numbers also support encoding and decoding specific IEEE 754+interchange formats via a 'Data.Binary.Binary' instance. -} module Numeric.Decimal        ( -- * Usage          -- $usage +         -- ** Advanced usage+         -- $advanced-usage+          -- * Arbitrary-precision decimal numbers          Decimal        , BasicDecimal        , ExtendedDecimal        , GeneralDecimal +         -- ** Number types with defined encodings+         -- $encodings+       , Decimal32+       , Decimal64+       , Decimal128+          -- ** Precision types-       , module Numeric.Decimal.Precision+       , Precision+       , FinitePrecision +       , P1 , P2 , P3 , P4 , P5 , P6 , P7 , P8 , P9 , P10+       , P11, P12, P13, P14, P15, P16, P17, P18, P19, P20+       , P21, P22, P23, P24, P25, P26, P27, P28, P29, P30+       , P31, P32, P33, P34, P35, P36, P37, P38, P39, P40+       , P41, P42, P43, P44, P45, P46, P47, P48, P49, P50++       , P75, P100, P150, P200, P250, P300, P400, P500, P1000, P2000++       , PPlus1, PTimes2++       , PInfinite++       , Pdecimal32+       , Pdecimal64+       , Pdecimal128+          -- ** Rounding types        , Rounding @@ -67,23 +96,14 @@          -- * Functions        , cast        , fromBool+       , fromOrdering        ) where +import Numeric.Decimal.Encoding import Numeric.Decimal.Number import Numeric.Decimal.Precision import Numeric.Decimal.Rounding --- | A decimal floating point number with 9 digits of precision, rounding half--- up-type BasicDecimal = Decimal P9 RoundHalfUp---- | A decimal floating point number with selectable precision, rounding half--- even-type ExtendedDecimal p = Decimal p RoundHalfEven---- | A decimal floating point number with infinite precision-type GeneralDecimal = ExtendedDecimal PInfinite- {- $usage  You should choose a decimal number type with appropriate precision and@@ -93,35 +113,52 @@ rounds half up.  * 'ExtendedDecimal' is a number type constructor with selectable precision-that rounds half even. For example, @'ExtendedDecimal' 'P34'@ is a number type-with 34 decimal digits of precision. There is a range of ready-made precisions-available, including 'P1' through 'P50' on up to 'P2000' (the IEEE 754-smallest and basic formats correspond to precisions 'P7', 'P16', or 'P34').-Alternatively, an arbitrary precision can be constructed through type-application of 'PPlus1' and/or 'PTimes2' to any existing precision.+that rounds half even. For example, @'ExtendedDecimal' 'P15'@ is a number type+with 15 decimal digits of precision. There is a range of ready-made precisions+available, including 'P1' through 'P50' on up to 'P2000'.  * 'GeneralDecimal' is a number type with infinite precision. Note that not all operations support numbers with infinite precision. +* 'Decimal32', 'Decimal64', and 'Decimal128' are specialized number types with+'Data.Binary.Binary' instances that implement the /decimal32/, /decimal64/,+and /decimal128/ interchange format encodings described in IEEE 754-2008.+These types have precisions of 7, 16, and 34 decimal digits, respectively, and+round half even.+ * The most versatile 'Decimal' type constructor is parameterized by both a precision and a rounding algorithm. For example, @'Decimal' 'P20' 'RoundDown'@ is a number type with 20 decimal digits of precision that rounds down (truncates). Several 'Rounding' algorithms are available to choose from. -It is suggested to create an alias for the type of numbers you wish to support-in your application. For example:--> type Number = ExtendedDecimal P16--A decimal number type may be used in a @default@ declaration, possibly-replacing 'Double' and/or 'Integer'. For example:+A decimal number type may be used in a @default@ declaration, for example+replacing 'Double':  > default (Integer, BasicDecimal)+-} -== Advanced usage+{- $advanced-usage -Additional operations and control beyond what is provided by the basic numeric-type classes are available through the use of "Numeric.Decimal.Arithmetic" and+Additional operations and control beyond what is provided by the standard type+classes are available through the use of "Numeric.Decimal.Arithmetic" and "Numeric.Decimal.Operation". Advanced string conversion is also available through "Numeric.Decimal.Conversion".++Arbitrary precisions can be constructed through type application of 'PPlus1'+and/or 'PTimes2' to any existing precision.++It is possible to create arbitrary width interchange format encodings with the+help of "Numeric.Decimal.Encoding".+-}++{- $encodings++These decimal number types have a 'Data.Binary.Binary' instance that+implements a specific interchange format encoding described in IEEE+754-2008. See "Numeric.Decimal.Encoding" for further details, including the+ability to create additional formats of arbitrary width.++Alternative rounding algorithms can be used through the more general 'Decimal'+type constructor and the special precision types 'Pdecimal32', 'Pdecimal64',+or 'Pdecimal128', e.g. @'Decimal' 'Pdecimal64' 'RoundCeiling'@. -}

@@ -71,8 +71,8 @@ -- arithmetic computation and manipulate its 'Context'.  -- | A context for decimal arithmetic, carrying signal flags, trap enabler--- state, and a trap handler, parameterized by precision @p@ and rounding--- algorithm @r@+-- state, and a trap handler, parameterized by precision @p@ and rounding mode+-- @r@ data Context p r =   Context { flags        :: Signals                             -- ^ The current signal flags of the context@@ -105,9 +105,8 @@         disabled = [Inexact, Rounded, Subnormal]  -- | Return a new context with all signal flags cleared, all traps disabled--- (IEEE 854 §7), using selectable precision (the IEEE 754 smallest and basic--- formats correspond to precisions 'P7', 'P16', or 'P34'), and rounding half--- even (IEEE 754 §4.3.3).+-- (IEEE 854 §7), using selectable precision, and rounding half even (IEEE 754+-- §4.3.3). extendedDefaultContext :: Context p RoundHalfEven extendedDefaultContext = newContext @@ -123,7 +122,7 @@   deriving Show  -- | A decimal arithmetic monad parameterized by the precision @p@ and--- rounding algorithm @r@+-- rounding mode @r@ newtype Arith p r a = Arith (ExceptT (Exception p r)                              (State (Context p r)) a) @@ -164,8 +163,8 @@ evalArith (Arith e) = evalState (runExceptT e)  -- | Perform a subcomputation using a different precision and/or rounding--- algorithm. The subcomputation is evaluated within a new context with all--- flags cleared and all traps disabled. Any flags set in the context of the+-- mode. The subcomputation is evaluated within a new context with all flags+-- cleared and all traps disabled. Any flags set in the context of the -- subcomputation are ignored, but if an exception is returned it will be -- re-raised within the current context. subArith :: Arith a b (Decimal a b) -> Arith p r (Decimal a b)@@ -181,7 +180,7 @@   where getPrecision' :: Precision p => p -> Arith p r (Maybe Int)         getPrecision' = return . precision --- | Return the rounding algorithm of the arithmetic context.+-- | Return the rounding mode of the arithmetic context. getRounding :: Rounding r => Arith p r RoundingAlgorithm getRounding = getRounding' undefined   where getRounding' :: Rounding r => r -> Arith p r RoundingAlgorithm

@@ -244,15 +244,18 @@               showString (replicate (fromIntegral $ -e - cl) '0') .               showString cs -  Inf  {             } -> showString "Infinity"-  QNaN { payload = p } -> showString  "NaN" . diag p-  SNaN { payload = p } -> showString "sNaN" . diag p+  Inf {                            } -> showString "Infinity"+  NaN { signaling = s, payload = p } -> sig s . showString  "NaN" . diag p    where signStr :: ShowS         signStr = case sign num of           Pos -> id           Neg -> showChar '-' +        sig :: Bool -> ShowS+        sig False = id+        sig True  = showChar 's'+         diag :: Payload -> ShowS         diag 0 = id         diag d = shows d@@ -295,34 +298,26 @@                         , exponent    = e - fromIntegral (length fracDigits)                         } -                digitsWithOptionalPoint = fractionalDigits <|> wholeDigits--                fractionalDigits = do-                  char '.'-                  fracDigits <- many1 parseDigit-                  return $ \e ->-                    Num { sign        = Pos-                        , coefficient = readDigits fracDigits-                        , exponent    = e - fromIntegral (length fracDigits)-                        }--                wholeDigits = do+                digitsWithOptionalPoint = do+                  fractional <- option False (char '.' *> pure True)                   digits <- many1 parseDigit+                  let offset | fractional = fromIntegral (length digits)+                             | otherwise  = 0                   return $ \e -> Num { sign        = Pos                                      , coefficient = readDigits digits-                                     , exponent    = e+                                     , exponent    = e - offset                                      }          parseExponentPart :: ReadP Exponent         parseExponentPart = do           parseString "E"-          parseSign negate <*> (readDigits <$> many1 parseDigit)+          parseSign negate <*> (fromIntegral . readDigits <$> many1 parseDigit)          parseInfinity :: ReadP (Decimal p r)         parseInfinity = do           parseString "Inf"           optional $ parseString "inity"-          return Inf { sign = Pos }+          return infinity          parseNaN :: ReadP (Decimal p r)         parseNaN = parseQNaN <|> parseSNaN@@ -330,13 +325,13 @@         parseQNaN :: ReadP (Decimal p r)         parseQNaN = do           p <- parseNaNPayload-          return QNaN { sign = Pos, payload = p }+          return qNaN { payload = p }          parseSNaN :: ReadP (Decimal p r)         parseSNaN = do           parseString "s"           p <- parseNaNPayload-          return SNaN { sign = Pos, payload = p }+          return sNaN { payload = p }          parseNaNPayload :: ReadP Payload         parseNaNPayload = do@@ -349,7 +344,7 @@         parseString :: String -> ReadP ()         parseString = mapM_ $ \c -> char (toLower c) <|> char (toUpper c) -        readDigits :: Num c => [Int] -> c+        readDigits :: [Int] -> Coefficient         readDigits = foldl' (\a b -> a * 10 + fromIntegral b) 0  {- $doctest-toNumber

@@ -0,0 +1,330 @@++{-# LANGUAGE FlexibleInstances #-}++-- | This module implements the decimal interchange format encodings described+-- in IEEE 754-2008, including the /decimal32/, /decimal64/, and /decimal128/+-- formats, as well as arbitrary width /decimal{k}/ formats through the use of+-- 'Format' with 'KPlus32' and\/or 'KTimes2'. For example, to use a+-- /decimal96/ format:+--+-- > type Decimal96 = ExtendedDecimal (Format (KPlus32 K64) DecimalCoefficient)+--+-- Currently only a decimal encoding of coefficients is implemented, but a+-- binary encoding may be added in the future.+module Numeric.Decimal.Encoding (+    -- * Primary convenience types+    Decimal32+  , Decimal64+  , Decimal128++    -- ** Precision types+  , Pdecimal32+  , Pdecimal64+  , Pdecimal128++    -- * Interchange format types+  , Format++    -- ** Format parameters+  , Parameters+  , K32+  , K64+  , K128+  , KPlus32+  , KTimes2++    -- ** Coefficient encodings+  -- , CoefficientEncoding+  , DecimalCoefficient+  , BinaryCoefficient+  ) where++import Prelude hiding (exponent)++import Data.Binary (Binary(get, put), Get)+import Data.Binary.Bits.Get (BitGet, getBool, getWord8, getWord16be, runBitGet)+import Data.Binary.Bits.Put (BitPut, putBool, putWord8, putWord16be, runBitPut)+import Data.Bits (bit, shiftL, shiftR, testBit, (.&.), (.|.))+import Data.Word (Word8, Word16)++import Numeric.Decimal.Number+import Numeric.Decimal.Precision++-- Decimal number types++-- | A decimal floating point number with 7 digits of precision, rounding half+-- even, and a 32-bit encoded representation using the /decimal32/ interchange+-- format (with a decimal encoding for the coefficient)+type Decimal32 = ExtendedDecimal Pdecimal32++-- | A decimal floating point number with 16 digits of precision, rounding+-- half even, and a 64-bit encoded representation using the /decimal64/+-- interchange format (with a decimal encoding for the coefficient)+type Decimal64 = ExtendedDecimal Pdecimal64++-- | A decimal floating point number with 34 digits of precision, rounding+-- half even, and a 128-bit encoded representation using the /decimal128/+-- interchange format (with a decimal encoding for the coefficient)+type Decimal128 = ExtendedDecimal Pdecimal128++-- Precision types++-- | A type with 'Precision' instance specifying /decimal32/ interchange+-- format parameters (using a decimal encoding for the coefficient) having an+-- effective precision of 7 decimal digits+type Pdecimal32 = Format K32 DecimalCoefficient++-- | A type with 'Precision' instance specifying /decimal64/ interchange+-- format parameters (using a decimal encoding for the coefficient) having an+-- effective precision of 16 decimal digits+type Pdecimal64 = Format K64 DecimalCoefficient++-- | A type with 'Precision' instance specifying /decimal128/ interchange+-- format parameters (using a decimal encoding for the coefficient) having an+-- effective precision of 34 decimal digits+type Pdecimal128 = Format K128 DecimalCoefficient++-- Format parameters++-- | Interchange format parameters used to define an encoding and derive the+-- format's /precision/ and E/max/+class Parameters k where+  -- | /k//32, the primary format parameter defining the encoding width as a+  -- multiple of 32 bits+  paramK32 :: k -> Int++-- | /p/, precision in digits+paramP :: Parameters k => k -> Int+paramP k = 9 * paramK32 k - 2++-- | /emax/+paramEmax :: Parameters k => k -> Exponent+paramEmax k = 3 * 2^(paramK32 k * 2 + 3)++-- | /bias/, /E/ − /q/+paramBias :: Parameters k => k -> Exponent+paramBias k = paramEmax k + fromIntegral (paramP k - 2)++-- | /w/, combination field width in bits − 5+paramW :: Parameters k => k -> Int+paramW k = paramK32 k * 2 + 4++-- | /t//10, trailing significand field width in 10-bit multiples+paramT10 :: Parameters k => k -> Int+paramT10 k = 3 * paramK32 k - 1++-- | Parameters for the /decimal32/ interchange format+data K32+instance Parameters K32 where+  paramK32 _ = 1++-- | Parameters for the /decimal64/ interchange format+type K64 = KPlus32 K32++-- | Parameters for the /decimal128/ interchange format+type K128 = KTimes2 K64++-- | Parameters for a /decimal{@k@ + 32}/ interchange format+data KPlus32 k+instance Parameters k => Parameters (KPlus32 k) where+  paramK32 t = paramK32 (minus32 t) + 1+    where minus32 :: KPlus32 k -> k+          minus32 = undefined++-- | Parameters for a /decimal{@k@ × 2}/ interchange format+data KTimes2 k+instance Parameters k => Parameters (KTimes2 k) where+  paramK32 t = paramK32 (div2 t) * 2+    where div2 :: KTimes2 k -> k+          div2 = undefined++-- | A class encapsulating coefficient encodings+class CoefficientEncoding c++-- | Specify a decimal encoding for the coefficient.+data DecimalCoefficient+instance CoefficientEncoding DecimalCoefficient++-- | Specify a binary encoding for the coefficient (currently unimplemented).+data BinaryCoefficient+instance CoefficientEncoding BinaryCoefficient++-- | A type (with a 'Precision' instance) for specifying interchange format+-- parameters @k@ and coefficient encoding @c@+data Format k c++formatK :: Format k c -> k+formatK = undefined++-- | This 'Precision' instance automatically computes the /precision/ and+-- E/max/ of decimal numbers that use this format.+instance Parameters k => Precision (Format k c) where+  precision = Just . paramP    . formatK+  eMax      = Just . paramEmax . formatK++instance Parameters k => FinitePrecision (Format k c)++-- | A 'Binary' instance is defined for interchange formats for which a+-- 'Parameters' instance exists, and covers particularly the 'Decimal32',+-- 'Decimal64', and 'Decimal128' types.+instance Parameters k => Binary (Decimal (Format k DecimalCoefficient) r) where+  put d = runBitPut $ putDecimal (paramW k) (paramT10 k) (paramBias k) d+    where k = formatK (decimalFormat d)++          decimalFormat :: Decimal (Format k c) r -> Format k c+          decimalFormat = undefined++  get = result+    where result = runBitGet $ getDecimal (paramW k) (paramT10 k) (paramBias k)+          k = formatK (getDecimalFormat result)++          getDecimalFormat :: Get (Decimal (Format k c) r) -> Format k c+          getDecimalFormat = undefined++-- Densely Packed Decimal++dpd2bcd :: Word16 -> (Word8, Word8, Word8)+dpd2bcd dpd = case mask 0 0xe of+  0xe -> case mask 4 0x6 of+    0x6 -> (    c7,     f4,     i0)+    0x4 -> (a9b8c7,     f4,     i0)+    0x2 -> (    c7, d9e8f4,     i0)+    _   -> (    c7,     f4, g9h8i0)+  0xc ->   (    c7, d6e5f4, g9h8i0)+  0xa ->   (a9b8c7,     f4, g6h5i0)+  0x8 ->   (a9b8c7, d6e5f4,     i0)+  _   ->   (a9b8c7, d6e5f4, g2h1i0)++  where a9b8c7 =              mask 7 7+        d6e5f4 =              mask 4 7+        d9e8f4 = mask 7 6 .|. mask 4 1+        g2h1i0 =              mask 0 7+        g6h5i0 = mask 4 6 .|. mask 0 1+        g9h8i0 = mask 7 6 .|. mask 0 1+        i0     = 8        .|. mask 0 1+        f4     = 8        .|. mask 4 1+        c7     = 8        .|. mask 7 1++        mask :: Int -> Word8 -> Word8+        mask s m = fromIntegral (shiftR dpd s) .&. m++bcd2dpd :: Word8 -> Word8 -> Word8 -> Word16+bcd2dpd d2 d1 d0 = case (d2 < 8, d1 < 8, d0 < 8) of+  (True , True , True ) ->      a9b8c7 .|.      d6e5f4          .|. g2h1i0+  (True , True , False) ->      a9b8c7 .|.      d6e5f4 .|. 0x08 .|.     i0+  (True , False, True ) ->      a9b8c7 .|. g6h5 .|. f4 .|. 0x0a .|.     i0+  (False, True , True ) -> g9h8 .|. c7 .|.      d6e5f4 .|. 0x0c .|.     i0+  (False, False, True ) -> g9h8 .|. c7 .|.          f4 .|. 0x0e .|.     i0+  (False, True , False) -> d9e8 .|. c7 .|.          f4 .|. 0x2e .|.     i0+  (True , False, False) ->      a9b8c7 .|.          f4 .|. 0x4e .|.     i0+  (False, False, False) ->          c7 .|.          f4 .|. 0x6e .|.     i0++  where a9b8c7 = isolate d2 7 7+        c7     = isolate d2 1 7+        d6e5f4 = isolate d1 7 4+        d9e8   = isolate d1 6 7+        f4     = isolate d1 1 4+        g2h1i0 = isolate d0 7 0+        g6h5   = isolate d0 6 4+        g9h8   = isolate d0 6 7+        i0     = isolate d0 1 0++        isolate :: Word8 -> Word8 -> Int -> Word16+        isolate d m = shiftL (fromIntegral $ d .&. m)++-- Low-level encoding/decoding++data CombinationField = Finite { exponentMSBs   :: Word8+                               , coefficientMSD :: Word8 }+                      | Infinity+                      | NotANumber++getCommon :: BitGet (Sign, CombinationField)+getCommon = do+  sign <- toEnum . fromEnum <$> getBool+  bits <- getWord8 5+  let cf = case bits of+        0x1e -> Infinity+        0x1f -> NotANumber+        _    -> let ab = shiftR bits 3 in case ab of+          0x03 -> Finite { exponentMSBs   = shiftR bits 1 .&. 0x03+                         , coefficientMSD = 0x08 .|. (bits .&. 0x01)+                         }+          _    -> Finite { exponentMSBs   = ab+                         , coefficientMSD = bits .&. 0x07+                         }+  return (sign, cf)++putCommon :: Sign -> CombinationField -> BitPut ()+putCommon sign cf = do+  putBool (toEnum . fromEnum $ sign)+  let bits = case cf of+        Finite { exponentMSBs = msbs, coefficientMSD = msd }+          | msd < 8   ->          shiftL msbs 3 .|.  msd+          | otherwise -> 0x18 .|. shiftL msbs 1 .|. (msd .&. 0x01)+        Infinity      -> 0x1e+        NotANumber    -> 0x1f+  putWord8 5 bits++getCoefficient :: CombinationField -> Int -> BitGet Coefficient+getCoefficient = getCoefficient' . getMSD++  where getCoefficient' :: Coefficient -> Int -> BitGet Coefficient+        getCoefficient' ic 0 = return ic+        getCoefficient' ic n = do+          (a, b, c) <- dpd2bcd <$> getWord16be 10+          let v = fromIntegral a * 100 + fromIntegral b * 10 + fromIntegral c+          getCoefficient' (ic * 1000 + v) (pred n)++        getMSD :: CombinationField -> Coefficient+        getMSD Finite { coefficientMSD = msd } = fromIntegral msd+        getMSD _                               = 0++getDecimal :: Int -> Int -> Exponent -> BitGet (Decimal p r)+getDecimal ecbits cclen bias = do+  (sign, cf) <- getCommon+  ec <- getWord16be ecbits+  coefficient <- getCoefficient cf cclen+  return $ case cf of+    Finite { exponentMSBs = msbs } ->+      let ee = shiftL (fromIntegral msbs) ecbits .|. fromIntegral ec+      in Num { sign = sign, coefficient = coefficient, exponent = ee - bias }+    Infinity -> Inf { sign = sign }+    NotANumber ->+      let s = testBit ec (ecbits - 1)+      in NaN { sign = sign, signaling = s, payload = coefficient }++putDecimal :: Int -> Int -> Exponent -> Decimal p r -> BitPut ()+putDecimal ecbits cclen bias x = do+  let msd : cc = digits x+      (cf, ee) = case x of+        Num { exponent = e } ->+          let cf = Finite { exponentMSBs   = fromIntegral (shiftR ee ecbits)+                          , coefficientMSD = msd+                          }+          in (cf, fromIntegral $ e + bias)+        Inf{} -> (Infinity, 0)+        NaN { signaling = s } -> (NotANumber, if s then bit (ecbits - 1) else 0)+  putCommon (sign x) cf+  putWord16be ecbits (ee .&. (bit ecbits - 1))+  putDigits cc++  where digits :: Decimal p r -> [Word8]+        digits x = let ds = case x of+                         Num { coefficient = c } -> digits' c+                         NaN { payload     = p } -> digits' p+                         Inf {                 } -> []+                   in replicate (1 + cclen * 3 - length ds) 0 ++ ds++        digits' :: Coefficient -> [Word8]+        digits' = go []+          where go ds 0 = ds+                go ds c = let (q, r) = c `quotRem` 10+                          in go (fromIntegral r : ds) q++        putDigits :: [Word8] -> BitPut ()+        putDigits (a : b : c : rest) = do+          putWord16be 10 (bcd2dpd a b c)+          putDigits rest+        putDigits [] = return ()+        putDigits _ = error "putDigits: invalid # digits"

@@ -0,0 +1,242 @@++module Numeric.Decimal.Exception (+  -- * Exceptional conditions+    clamped+  , conversionSyntax+  , divisionByZero+  , divisionImpossible+  , divisionUndefined+  , inexact+  , insufficientStorage+  , invalidContext+  , invalidOperation+  , overflow+  , rounded+  , subnormal+  , underflow+  ) where++import Prelude hiding (exponent)++import Numeric.Decimal.Arithmetic+import Numeric.Decimal.Number+import Numeric.Decimal.Precision+import Numeric.Decimal.Rounding++-- | This occurs and signals 'Clamped' if the exponent of a result has been+-- altered in order to fit the constraints of a specific concrete+-- representation. This may occur when the exponent of a zero result would be+-- outside the bounds of a representation, or (in the IEEE 754 interchange+-- formats) when a large normal number would have an encoded exponent that+-- cannot be represented. In this latter case, the exponent is reduced to fit+-- and the corresponding number of zero digits are appended to the coefficient+-- (“fold-down”). The condition always occurs when a subnormal value rounds to+-- zero.+clamped :: Decimal p r -> Arith p r (Decimal p r)+clamped = raiseSignal Clamped++-- | This occurs and signals 'InvalidOperation' if a string is being converted+-- to a number and it does not conform to the numeric string syntax. The+-- result is @[0,qNaN]@.+conversionSyntax :: Arith p r (Decimal p r)+conversionSyntax = raiseSignal InvalidOperation qNaN++-- | This occurs and signals 'DivisionByZero' if division of a finite number+-- by zero was attempted (during a 'Numeric.Decimal.Operation.divideInteger'+-- or 'Numeric.Decimal.Operation.divide' operation, or a+-- 'Numeric.Decimal.Operation.power' operation with negative right-hand+-- operand), and the dividend was not zero.+--+-- The result of the operation is @[@/sign/@,inf]@, where /sign/ is the+-- exclusive or of the signs of the operands for divide, or is 1 for an odd+-- power of −0, for power.+divisionByZero :: Decimal p r -> Arith p r (Decimal p r)+divisionByZero = raiseSignal DivisionByZero++-- | This occurs and signals 'InvalidOperation' if the integer result of a+-- 'Numeric.Decimal.Operation.divideInteger' or+-- 'Numeric.Decimal.Operation.remainder' operation had too many digits (would+-- be longer than /precision/). The result is @[0,qNaN]@.+divisionImpossible :: Arith p r (Decimal p r)+divisionImpossible = raiseSignal InvalidOperation qNaN++-- | This occurs and signals 'InvalidOperation' if division by zero was+-- attempted (during a 'Numeric.Decimal.Operation.divideInteger',+-- 'Numeric.Decimal.Operation.divide', or+-- 'Numeric.Decimal.Operation.remainder' operation), and the dividend is also+-- zero. The result is @[0,qNaN]@.+divisionUndefined :: Arith p r (Decimal p r)+divisionUndefined = raiseSignal InvalidOperation qNaN++-- | This occurs and signals 'Inexact' whenever the result of an operation is+-- not exact (that is, it needed to be rounded and any discarded digits were+-- non-zero), or if an overflow or underflow condition occurs. The result in+-- all cases is unchanged.+--+-- The 'Inexact' signal may be tested (or trapped) to determine if a given+-- operation (or sequence of operations) was inexact.+inexact :: Decimal p r -> Arith p r (Decimal p r)+inexact = raiseSignal Inexact++-- | For many implementations, storage is needed for calculations and+-- intermediate results, and on occasion an arithmetic operation may fail due+-- to lack of storage. This is considered an operating environment error,+-- which can be either be handled as appropriate for the environment, or+-- treated as an Invalid operation condition. The result is @[0,qNaN]@.+insufficientStorage :: Arith p r (Decimal p r)+insufficientStorage = invalidOperation qNaN++-- | This occurs and signals 'InvalidOperation' if an invalid context was+-- detected during an operation. This can occur if contexts are not checked on+-- creation and either the /precision/ exceeds the capability of the+-- underlying concrete representation or an unknown or unsupported /rounding/+-- was specified. These aspects of the context need only be checked when the+-- values are required to be used. The result is @[0,qNaN]@.+invalidContext :: Arith p r (Decimal p r)+invalidContext = raiseSignal InvalidOperation qNaN++-- | This occurs and signals 'InvalidOperation' if:+--+-- * an operand to an operation is @[s,sNaN]@ or @[s,sNaN,d]@ (any /signaling/+-- NaN)+--+-- * an attempt is made to add @[0,inf]@ to @[1,inf]@ during an addition or+-- subtraction operation+--+-- * an attempt is made to multiply 0 by @[0,inf]@ or @[1,inf]@+--+-- * an attempt is made to divide either @[0,inf]@ or @[1,inf]@ by either+-- @[0,inf]@ or @[1,inf]@+--+-- * the divisor for a remainder operation is zero+--+-- * the dividend for a remainder operation is either @[0,inf]@ or @[1,inf]@+--+-- * either operand of the 'Numeric.Decimal.Operation.quantize' operation is+-- infinite, or the result of a 'Numeric.Decimal.Operation.quantize' operation+-- would require greater precision than is available+--+-- * the operand of the 'Numeric.Decimal.Operation.ln' or the+-- 'Numeric.Decimal.Operation.log10' operation is less than zero+--+-- * the operand of the 'Numeric.Decimal.Operation.squareRoot' operation has a+-- /sign/ of 1 and a non-zero /coefficient/+--+-- * both operands of the 'Numeric.Decimal.Operation.power' operation are+-- zero, or if the left-hand operand is less than zero and the right-hand+-- operand does not have an integral value or is infinite+--+-- * an operand is invalid; for example, certain values of concrete+-- representations may not correspond to numbers — an implementation is+-- permitted (but is not required) to detect these invalid values and raise+-- this condition.+--+-- The result of the operation after any of these invalid operations is+-- @[0,qNaN]@ except when the cause is a signaling NaN, in which case the+-- result is @[s,qNaN]@ or @[s,qNaN,d]@ where the sign and diagnostic are+-- copied from the signaling NaN.+invalidOperation :: Decimal a b -> Arith p r (Decimal p r)+invalidOperation n = raiseSignal InvalidOperation $ case n of+  NaN { signaling = True } -> n { signaling = False }+  _                        -> qNaN++-- | This occurs and signals 'Overflow' if the /adjusted exponent/ of a result+-- (from a conversion or from an operation that is not an attempt to divide by+-- zero), after rounding, would be greater than the largest value that can be+-- handled by the implementation (the value E/max/).+--+-- The result depends on the rounding mode:+--+-- * For 'RoundHalfUp' and 'RoundHalfEven' (and for 'RoundHalfDown' and+-- 'RoundUp', if implemented), the result of the operation is [sign,@inf@],+-- where /sign/ is the sign of the intermediate result.+--+-- * For 'RoundDown', (and 'Round05Up', if implemented), the result is the+-- largest finite number that can be represented in the current /precision/,+-- with the sign of the intermediate result.+--+-- * For 'RoundCeiling', the result is the same as for 'RoundDown' if the sign+-- of the intermediate result is 1, or is @[0,inf]@ otherwise.+--+-- * For 'RoundFloor', the result is the same as for 'RoundDown' if the sign+-- of the intermediate result is 0, or is @[1,inf]@ otherwise.+--+-- In all cases, 'inexact' and 'rounded' will also be raised.+--+-- Note: IEEE 854 §7.3 requires that the result delivered to a trap handler be+-- different, depending on whether the overflow was the result of a conversion+-- or of an arithmetic operation. This specification deviates from IEEE 854 in+-- this respect; however, an implementation could comply with IEEE 854 by+-- providing a separate mechanism for the special result to a trap+-- handler. IEEE 754 has no such requirement.+overflow :: (Precision p, Rounding r) => Decimal p r -> Arith p r (Decimal p r)+overflow ir = result >>= raiseSignal Overflow >>= inexact >>= rounded++  where result :: (Precision p, Rounding r) => Arith p r (Decimal p r)+        result = getRounding >>= \r -> case r of+          RoundHalfUp   -> signedInfinity+          RoundHalfEven -> signedInfinity+          RoundHalfDown -> signedInfinity+          RoundUp       -> signedInfinity++          RoundDown     -> largestFinite+          Round05Up     -> largestFinite++          RoundCeiling  -> case sign ir of+            Neg -> largestFinite+            Pos -> return infinity++          RoundFloor    -> case sign ir of+            Pos -> largestFinite+            Neg -> return infinity { sign = Neg }++        signedInfinity :: Arith p r (Decimal p r)+        signedInfinity = return infinity { sign = sign ir }++        largestFinite :: Precision p => Arith p r (Decimal p r)+        largestFinite = getPrecision >>= \p ->+          let x = Num { sign = sign ir+                      , coefficient = undefined -- 10^p - 1+                      , exponent = undefined -- eMax x+                      }+          in return x++-- | This occurs and signals 'Rounded' whenever the result of an operation is+-- rounded (that is, some zero or non-zero digits were discarded from the+-- coefficient), or if an overflow or underflow condition occurs. The result+-- in all cases is unchanged.+--+-- The 'Rounded' signal may be tested (or trapped) to determine if a given+-- operation (or sequence of operations) caused a loss of precision.+rounded :: Decimal p r -> Arith p r (Decimal p r)+rounded = raiseSignal Rounded++-- | This occurs and signals 'Subnormal' whenever the result of a conversion+-- or operation is subnormal (that is, its adjusted exponent is less than+-- E/min/, before any rounding). The result in all cases is unchanged.+--+-- The 'Subnormal' signal may be tested (or trapped) to determine if a given+-- or operation (or sequence of operations) yielded a subnormal result.+subnormal :: Decimal p r -> Arith p r (Decimal p r)+subnormal = raiseSignal Subnormal++-- | This occurs and signals 'Underflow' if a result is inexact and the+-- /adjusted exponent/ of the result would be smaller (more negative) than the+-- smallest value that can be handled by the implementation (the value+-- E/min/). That is, the result is both inexact and subnormal.+--+-- The result after an underflow will be a subnormal number rounded, if+-- necessary, so that its exponent is not less than E/tiny/. This may result+-- in 0 with the sign of the intermediate result and an exponent of E/tiny/.+--+-- In all cases, 'inexact', 'rounded', and 'subnormal' will also be raised.+--+-- Note: IEEE 854 §7.4 requires that the result delivered to a trap handler be+-- different, depending on whether the underflow was the result of a+-- conversion or of an arithmetic operation. This specification deviates from+-- IEEE 854 in this respect; however, an implementation could comply with IEEE+-- 854 by providing a separate mechanism for the result to a trap+-- handler. IEEE 754 has no such requirement.+underflow :: Decimal p r -> Arith p r (Decimal p r)+underflow x = undefined >>= raiseSignal Underflow >>=+  subnormal >>= inexact >>= rounded

@@ -0,0 +1,11 @@++module Numeric.Decimal.Exception (+    inexact+  , rounded+  ) where++import {-# SOURCE #-} Numeric.Decimal.Arithmetic+import {-# SOURCE #-} Numeric.Decimal.Number++inexact :: Decimal p r -> Arith p r (Decimal p r)+rounded :: Decimal p r -> Arith p r (Decimal p r)

@@ -11,6 +11,10 @@        , Payload         , Decimal(..)+       , BasicDecimal+       , ExtendedDecimal+       , GeneralDecimal+        , zero        , oneHalf        , one@@ -23,7 +27,9 @@         , flipSign        , cast+        , fromBool+       , fromOrdering         , isPositive        , isNegative@@ -38,6 +44,7 @@  import Prelude hiding (exponent) +import Control.DeepSeq (NFData(..)) import Control.Monad (join) import Data.Bits (Bits(..), FiniteBits(..)) import Data.Char (isSpace)@@ -55,14 +62,13 @@  import qualified GHC.Real -{- $setup->>> :load Harness--}- data Sign = Pos  -- ^ Positive or non-negative           | Neg  -- ^ Negative           deriving (Eq, Enum) +instance NFData Sign where+  rnf s = s `seq` ()+ negateSign :: Sign -> Sign negateSign Pos = Neg negateSign Neg = Pos@@ -83,25 +89,34 @@   _  -> Pos  type Coefficient = Natural-type Exponent    = Int+type Exponent    = Integer type Payload     = Coefficient  -- | A decimal floating point number with selectable precision and rounding -- algorithm data Decimal p r-  = Num  { sign        :: Sign-         , coefficient :: Coefficient-         , exponent    :: Exponent-         }-  | Inf  { sign        :: Sign-         }-  | QNaN { sign        :: Sign-         , payload     :: Payload-         }-  | SNaN { sign        :: Sign-         , payload     :: Payload-         }+  = Num { sign        :: Sign+        , coefficient :: Coefficient+        , exponent    :: Exponent+        }+  | Inf { sign        :: Sign+        }+  | NaN { sign        :: Sign+        , signaling   :: Bool+        , payload     :: Payload+        } +-- | A decimal floating point number with 9 digits of precision, rounding half+-- up+type BasicDecimal = Decimal P9 RoundHalfUp++-- | A decimal floating point number with selectable precision, rounding half+-- even+type ExtendedDecimal p = Decimal p RoundHalfEven++-- | A decimal floating point number with infinite precision+type GeneralDecimal = ExtendedDecimal PInfinite+ -- | The 'Show' instance uses the 'toScientificString' operation from -- "Numeric.Decimal.Conversion". instance Show (Decimal p r) where@@ -115,67 +130,61 @@                     | (n, s) <- readParen False                       (readP_to_S toNumber . dropWhile isSpace) str ] -{- $doctest-Read->>> fmap toRep (read "Just 123" :: Maybe GeneralDecimal)-Just (N (0,123,0))-->>> fmap toRep (read "Just (-12.0)" :: Maybe GeneralDecimal)-Just (N (1,120,-1))--}+decimalPrecision :: Decimal p r -> p+decimalPrecision = undefined  instance Precision p => Precision (Decimal p r) where   precision = precision . decimalPrecision-    where decimalPrecision :: Decimal p r -> p-          decimalPrecision = undefined+  eMax      = eMax      . decimalPrecision+  eMin      = eMin      . decimalPrecision -evalOp :: Arith p r (Decimal p r) -> Decimal p r-evalOp op = either exceptionResult id $ evalArith op newContext+-- This assumes the arithmetic operation does not trap any signals, which+-- could result in an exception being thrown (and returned in a Left value).+evalOp :: Arith p r a -> a+evalOp op = let Right r = evalArith op newContext in r -type GeneralDecimal = Decimal PInfinite RoundHalfEven+evalOp' :: Arith p RoundHalfEven a -> a+evalOp' = evalOp +compareDecimal :: Decimal a b -> Decimal c d -> Either GeneralDecimal Ordering+compareDecimal x y = evalOp (x `Op.compare` y)++-- | Note that NaN values are not equal to any value, including other NaNs. instance Eq (Decimal p r) where-  x == y = case evalOp (x `Op.compare` y) :: GeneralDecimal of-    Num { coefficient = 0 } -> True-    _                       -> False+  x == y = case compareDecimal x y of+    Right EQ -> True+    _        -> False -instance (Precision p, Rounding r) => Ord (Decimal p r) where-  x `compare` y = case evalOp (x `Op.compare` y) :: GeneralDecimal of-    Num { coefficient = 0 } -> EQ-    Num { sign = Neg      } -> LT-    Num { sign = Pos      } -> GT-    _                       -> GT  -- match Prelude behavior for NaN+-- | Unlike the instances for 'Float' and 'Double', the 'compare' method in+-- this instance uses a total ordering over all possible values. Note that+-- @'compare' x y == 'EQ'@ does not imply @x == y@ (and similarly for 'LT' and+-- 'GT') in the cases where @x@ or @y@ are NaN values.+instance Ord (Decimal p r) where+  compare x y = case compareDecimal x y of+    Right o -> o+    Left _  -> evalOp (x `Op.compareTotal` y) -  x < y = case evalOp (x `Op.compare` y) :: GeneralDecimal of-    Num { sign = Neg      } -> True-    _                       -> False+  x < y = case compareDecimal x y of+    Right LT -> True+    _        -> False -  x <= y = case evalOp (x `Op.compare` y) :: GeneralDecimal of-    Num { sign = Neg      } -> True-    Num { coefficient = 0 } -> True-    _                       -> False+  x <= y = case compareDecimal x y of+    Right LT -> True+    Right EQ -> True+    _        -> False -  x > y = case evalOp (x `Op.compare` y) :: GeneralDecimal of-    Num { coefficient = 0 } -> False-    Num { sign = Pos      } -> True-    _                       -> False+  x > y = case compareDecimal x y of+    Right GT -> True+    _        -> False -  x >= y = case evalOp (x `Op.compare` y) :: GeneralDecimal of-    Num { sign = Pos      } -> True-    _                       -> False+  x >= y = case compareDecimal x y of+    Right GT -> True+    Right EQ -> True+    _        -> False    max x y = evalOp (Op.max x y)   min x y = evalOp (Op.min x y) -{- $doctest-Ord-prop> x > y ==> max x y == x && max y x == (x :: BasicDecimal)-prop> x < y ==> min x y == x && min y x == (x :: BasicDecimal)--prop> max x y == x ==> x >= y-prop> max x y == y ==> y >= x-prop> min x y == x ==> x <= y-prop> min x y == y ==> y <= x--}- -- | Unlike the instances for 'Float' and 'Double', the lists returned by the -- 'enumFromTo' and 'enumFromThenTo' methods in this instance terminate with -- the last element strictly less than (greater than in the case of a negative@@ -204,17 +213,6 @@              => Decimal p r -> Decimal p r -> [Decimal p r] enumFromWith x i = x : enumFromWith (x + i) i -{- $doctest-Enum->>> [0, 0.1 .. 2] :: [BasicDecimal]-[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0]-->>> [2, 1.9 .. 0] :: [BasicDecimal]-[2,1.9,1.8,1.7,1.6,1.5,1.4,1.3,1.2,1.1,1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1,0.0]-->>> [1.7 .. 5.7] :: [BasicDecimal]-[1.7,2.7,3.7,4.7,5.7]--}- instance (Precision p, Rounding r) => Num (Decimal p r) where   x + y = evalOp (x `Op.add`      y)   x - y = evalOp (x `Op.subtract` y)@@ -224,10 +222,11 @@   abs    = evalOp . Op.abs    signum n = case n of-    Num { coefficient = 0 } -> zero-    Num { sign = s        } -> one { sign = s }-    Inf { sign = s        } -> one { sign = s }-    _                       -> n+    Num { sign = s, coefficient = c }+      | c == 0       -> zero { sign = s }+      | otherwise    -> one  { sign = s }+    Inf { sign = s } -> one  { sign = s }+    _                -> qNaN    fromInteger x = cast     Num { sign        = signMatch x@@ -235,15 +234,6 @@         , exponent    = 0         } -{- $doctest-Num-prop> x + x == x * (2 :: GeneralDecimal)-prop> isFinite x ==> x - x == (0 :: GeneralDecimal)-prop> isFinite x ==> x + negate x == (0 :: GeneralDecimal)-prop> abs x >= (0 :: GeneralDecimal)--prop> abs x * signum x == (x :: GeneralDecimal)--}- instance (Precision p, Rounding r) => Real (Decimal p r) where   toRational Num { sign = s, coefficient = c, exponent = e }     | e >= 0    = fromInteger $ signFunc s (fromIntegral c * 10^e)@@ -258,16 +248,6 @@                        d = fromInteger (denominator r) :: GeneralDecimal                    in evalOp (n `Op.divide` d) -{- $doctest-Fractional-prop> (4.14 :: Decimal P2 RoundHalfUp)   == 4.1-prop> (4.15 :: Decimal P2 RoundHalfUp)   == 4.2-prop> (4.15 :: Decimal P2 RoundHalfDown) == 4.1-prop> (4.15 :: Decimal P2 RoundHalfEven) == 4.2-prop> (4.25 :: Decimal P2 RoundHalfEven) == 4.2-prop> (4.35 :: Decimal P2 RoundHalfEven) == 4.4-prop> (4.45 :: Decimal P2 RoundHalfEven) == 4.4--}- instance (FinitePrecision p, Rounding r) => RealFrac (Decimal p r) where   properFraction x@Num { sign = s, coefficient = c, exponent = e }     | e < 0     = (n, f)@@ -277,95 +257,254 @@           (q, r) = c `quotRem` (10^(-e))   properFraction nan = (0, nan) -{- $doctest-RealFrac-prop> let (n,f) = properFraction (x :: BasicDecimal) in x == fromIntegral n + f-prop> let (n,f) = properFraction (x :: BasicDecimal) in (x < 0 && n <= 0) || (x >= 0 && n >= 0)-prop> let (n,f) = properFraction (x :: BasicDecimal) in (x < 0 && f <= 0) || (x >= 0 && f >= 0)-prop> let (n,f) = properFraction (x :: BasicDecimal) in isFinite f ==> abs f < 1--}+-- | Compute a generalized continued fraction to maximum precision. A hint is+-- used to indicate the minimum number of terms that should be generated+-- before (expensively) examining the results for convergence.+continuedFraction :: FinitePrecision p+                  => Int -> Integer -> [(Integer, Integer)] -> ExtendedDecimal p+continuedFraction m b0 ((a1, b1) : ps) = convergent (max 0 $ m - 2) x0 x1 x1' ps+  where x0  = (aa0, bb0)+        x1  = (aa1, bb1)+        x1' = fromRational (aa1 % bb1)+        aa0 = b0+        bb0 = 1+        aa1 = b1 * b0 + a1+        bb1 = b1 +        convergent m (aa0, bb0) x1@(aa1, bb1) x1' ((a2, b2) : ps)+          | m == 0 && x2' == x1' = x2'+          | otherwise            = convergent (max 0 $ m - 1) x1 x2 x2' ps+          where x2  = (aa2, bb2)+                x2' = fromRational (aa2 % bb2)+                aa2 = b2 * aa1 + a2 * aa0+                bb2 = b2 * bb1 + a2 * bb0+        convergent _ _ _ x [] = x++continuedFraction _ b0 [] = fromInteger b0+ -- | Compute an infinite series to maximum precision.-infiniteSeries :: (FinitePrecision p, Rounding r)-               => [Decimal p r] -> Decimal p r-infiniteSeries = series zero+infiniteSeries :: FinitePrecision p+               => (ExtendedDecimal p -> ExtendedDecimal p -> ExtendedDecimal p)+               -> [ExtendedDecimal p] -> ExtendedDecimal p+infiniteSeries op ~(x:xs) = series x xs   where series n (x:xs)           | n' == n   = n'           | otherwise = series n' xs-          where n' = n + x+          where n' = n `op` x         series n []   = n --- | Compute the arcsine of the argument to maximum precision using series--- expansion.-arcsine :: (FinitePrecision p, Rounding r) => Decimal p r -> Decimal p r-arcsine x = infiniteSeries (x : series 1 2 x 3)-  where series n d x i =-          let x' = x * x2-          in (n * x') / (d * i) : series (n * i) (d * (i + one)) x' (i + two)+-- | Compute the inverse tangent of the argument to maximum precision using+-- series expansion.+seriesArctan :: FinitePrecision p => Decimal p r -> ExtendedDecimal p+seriesArctan z = infiniteSeries (+) (z' : series True three z')+  where series neg d z =+          let z' = z * z2+              n | neg       = flipSign z'+                | otherwise = z'+          in (n / d) : series (not neg) (d + two) z'+        z' = castRounding z+        z2 = z' * z'++-- | Compute the inverse sine of the argument to maximum precision using+-- series expansion.+seriesArcsin :: FinitePrecision p => Decimal p r -> ExtendedDecimal p+seriesArcsin z = infiniteSeries (+) (z' : series one two z' three)+  where series n d z i =+          let z' = z * z2+          in (n * z') / (d * i) : series (n * i) (d * (i + one)) z' (i + two)+        z' = castRounding z+        z2 = z' * z'++-- | Compute the inverse tangent of the (decimal) argument to maximum+-- precision.+arctan :: (FinitePrecision p, Rounding r) => Decimal p r -> ExtendedDecimal p+arctan z@Num {          } = arctan' (toRational z)+arctan   Inf { sign = s } = signFunc s halfPi+arctan   _                = qNaN++-- | Compute the inverse tangent of the (rational) argument to maximum+-- precision using a generalized continued fraction.+arctan' :: FinitePrecision p => Rational -> ExtendedDecimal p+arctan' z = continuedFraction m 0 $ (x, y) : partials 1 0 y+  where x = numerator   z+        y = denominator z+        m = fromInteger $ (42 * abs x) `div` y  -- estimated minimum # terms+         x2 = x * x+        ty = 2 * y --- | Compute π to maximum precision using the arcsine series expansion.-seriesPi :: FinitePrecision p => Decimal p RoundHalfEven-seriesPi = 6 * arcsine oneHalf+        -- [ (nx * nx, (n * 2 + 1) * y) | n <- [1..], let nx = n * x ]+        partials n a b =+          let a' = a + n * x2+              b' = b +     ty+          in (a', b') : partials (n + 2) a' b' +-- | Compute the inverse sine of the argument to maximum precision.+arcsin :: FinitePrecision p => Decimal p r -> ExtendedDecimal p+arcsin z = let z' = castUp z+           in castDown' $ two * arctan (z' / (one + sqrt (one - z' * z')))++-- | Compute the inverse cosine of the argument to maximum precision.+arccos :: FinitePrecision p => Decimal p r -> ExtendedDecimal p+arccos = castDown' . (halfPi -) . arcsin . castUp++-- | Compute π to maximum precision using the inverse sine series expansion.+seriesPi :: FinitePrecision p => ExtendedDecimal p+seriesPi = castDown' $ 6 * arcsin oneHalf++-- | Compute π to maximum precision using the best-known Machin-like formula.+machinPi :: FinitePrecision p => ExtendedDecimal p+machinPi = castDown' $ 16 * arctan' (1 % 5) - 4 * arctan' (1 % 239)++-- | Compute π to maximum precision using a generalized continued fraction+-- that converges linearly, adding at least three decimal digits of precision+-- per four terms.+cfPi :: FinitePrecision p => ExtendedDecimal p+cfPi = pi'+  where pi'    = continuedFraction m+                 0 $ (4, 1) : [ (n * n, n * 2 + 1) | n <- [1..] ]+        Just p = precision pi'+        m      = (p `div` 3) * 4++-- | Precomputed π to a precision of 50 digits+fastPi :: FinitePrecision p => ExtendedDecimal p+fastPi = 3.1415926535897932384626433832795028841971693993751++-- | Compute π/2 to maximum precision.+halfPi :: FinitePrecision p => ExtendedDecimal p+halfPi = castDown' $ pi * oneHalf++-- | Compute π/4 to maximum precision.+quarterPi :: FinitePrecision p => ExtendedDecimal p+quarterPi = castDown' $ pi * oneQuarter++-- | Compute (cos 𝛽, sin 𝛽) to maximum precision using Volder's algorithm+-- (CORDIC).+cordic :: FinitePrecision p+       => ExtendedDecimal p -> (ExtendedDecimal p, ExtendedDecimal p)+cordic beta@Num{}+  | beta >          halfPi = negatePair $ cordic (beta - pi)+  | beta < flipSign halfPi = negatePair $ cordic (beta + pi)+  | isZero beta            = (one, zero)+  | otherwise              = cordic' beta (one, zero) one angles++  where negatePair (x, y) = (flipSign x, flipSign y)++        angles = quarterPi : [ arctan' z | let half = 1 % 2+                                         , z <- iterate (* half) half ]++        cordic' beta v@(x, y) powerOfTwo ~(angle:angles)+          | v' == v   = (k * x, k * y)+          | otherwise = cordic' beta' v' powerOfTwo' angles+          where isNegBeta = isNegative beta+                beta'  | isNegBeta = beta + angle+                       | otherwise = beta - angle+                factor | isNegBeta = powerOfTwo { sign = Neg }+                       | otherwise = powerOfTwo+                v' = (x - factor * y, factor * x + y)+                powerOfTwo' = powerOfTwo * oneHalf++        -- K = lim {n→∞} K(n)+        -- K(n) = prod {i=0..n-1} 1 / sqrt (1 + 2^(-2 * i))+        k | p <= 50   = fastK+          | otherwise = seriesK+          where Just p  = precision k+                fastK   = 0.60725293500888125616944675250492826311239085215009+                seriesK = infiniteSeries (*)+                  [ recip $ sqrt (one + x) | x <- iterate (* oneQuarter) one ]++cordic _ = (qNaN, qNaN)++-- | Cast a number to a number with two additional digits of precision and+-- rounding half even.+castUp :: Precision p => Decimal p r -> ExtendedDecimal (PPlus1 (PPlus1 p))+castUp = coerce+ -- | Cast a number with two additional digits of precision down to a number--- with the desired precision.-castDown :: (Precision p, Rounding r)-         => Decimal (PPlus1 (PPlus1 p)) a -> Decimal p r-castDown = cast+-- with the desired precision, rounding half even.+castDown' :: Precision p+          => ExtendedDecimal (PPlus1 (PPlus1 p)) -> ExtendedDecimal p+castDown' = cast -notyet :: String -> a-notyet = error . (++ ": not yet implemented")+-- | Cast a number with two additional digits of precision down to a number+-- with the desired precision, rounding half even, but returning a number type+-- with arbitrary rounding.+castDown :: (Precision p, Rounding r)+         => ExtendedDecimal (PPlus1 (PPlus1 p)) -> Decimal p r+castDown = castRounding . castDown' --- | The trigonometric and hyperbolic 'Floating' methods (other than the--- precision-dependent constant 'pi') are not yet implemented.+-- | The constant 'pi' is precision-dependent. instance (FinitePrecision p, Rounding r) => Floating (Decimal p r) where-  pi = castDown seriesPi+  pi = castRounding pi'+    where pi' | p <= 50   = fastPi+              | otherwise = cfPi+          Just p = precision pi'    exp = castRounding . evalOp . Op.exp   log = castRounding . evalOp . Op.ln -  logBase 10 x = castRounding $ evalOp (Op.log10 x)-  logBase _  1 = zero-  logBase b  x = evalOp (join $ Op.divide <$> Op.ln x <*> Op.ln b)+  logBase b@Num{} x+    | b == ten = castRounding $ evalOp (Op.log10 x)+    | x == one = case b `compare` one of+        LT -> zero { sign = Neg }+        EQ -> qNaN+        GT -> zero+    | x == b && not (isZero b) = one+  logBase b x = evalOp (join $ Op.divide <$> Op.ln x <*> Op.ln b)    x ** y = evalOp (x `Op.power` y)    sqrt = castRounding . evalOp . Op.squareRoot -  sin   = notyet "sin"-  cos   = notyet "cos"--  asin  = notyet "asin"-  acos  = notyet "acos"-  atan  = notyet "atan"--  sinh  = notyet "sinh"-  cosh  = notyet "cosh"--  asinh = notyet "asinh"-  acosh = notyet "acosh"-  atanh = notyet "atanh"+  sin  = castDown . snd                . cordic . castUp+  cos  = castDown . fst                . cordic . castUp+  tan  = castDown . uncurry (flip (/)) . cordic . castUp -{- $doctest-Floating-prop> realToFrac (pi :: ExtendedDecimal P16) == (pi :: Double)+  asin = castRounding . arcsin+  acos = castRounding . arccos+  atan = castRounding . arctan -prop> y >= 0 ==> (x :: BasicDecimal) ** fromInteger y == x ^ y+  -- sinh x = let ex = exp x in (ex^2 - 1) / (2 * ex)+  sinh x = castDown . evalOp' $+    Op.exp x >>= \ex -> two `Op.multiply` ex >>= \tex ->+    ex `Op.multiply` ex >>= (`Op.subtract` one) >>= (`Op.divide` tex)+  -- cosh x = let ex = exp x in (ex^2 + 1) / (2 * ex)+  cosh x = castDown . evalOp' $+    Op.exp x >>= \ex -> two `Op.multiply` ex >>= \tex ->+    ex `Op.multiply` ex >>= (`Op.add` one) >>= (`Op.divide` tex)+  -- tanh x = let e2x = exp (2 * x) in (e2x - 1) / (e2x + 1)+  tanh x = castDown . evalOp' $+    two `Op.multiply` x >>= Op.exp >>= \e2x ->+    e2x `Op.subtract` one >>= \e2xm1 -> e2x `Op.add` one >>= (e2xm1 `Op.divide`) -prop> isFinite x && x >= 0 ==> coefficient (sqrt (x * x) - (x :: ExtendedDecimal P16)) <= 1--}+  -- asinh x = log (x + sqrt (x^2 + 1))+  asinh x = castDown . evalOp' $ x `Op.multiply` x >>=+    (`Op.add` one) >>= Op.squareRoot >>= (x `Op.add`) >>= Op.ln+  -- acosh x = log (x + sqrt (x^2 - 1))+  acosh x = castDown . evalOp' $ x `Op.multiply` x >>=+    (`Op.subtract` one) >>= Op.squareRoot >>= (x `Op.add`) >>= Op.ln+  -- atanh x = log ((1 + x) / (1 - x)) / 2+  atanh x = castDown . evalOp' $ one `Op.add` x >>= \xp1 ->+    one `Op.subtract` x >>= (xp1 `Op.divide`) >>= Op.ln >>= Op.multiply oneHalf  instance (FinitePrecision p, Rounding r) => RealFloat (Decimal p r) where   floatRadix  _ = 10   floatDigits x = let Just p = precision x in p-  floatRange  _ = (minBound, maxBound)  -- ?+  floatRange  x = let Just emin = eMin x+                      Just emax = eMax x+                  in (fromIntegral emin, fromIntegral emax)    decodeFloat x = case x of-    Num  { sign = s, coefficient = c, exponent = e } -> (m, n)-      where m = signFunc s (fromIntegral c)-            n = fromIntegral e-    Inf  { sign = s              } -> (special s 0, maxBound    )-    QNaN { sign = s, payload = p } -> (special s p, minBound    )-    SNaN { sign = s, payload = p } -> (special s p, minBound + 1)+    Num { sign = s, coefficient = c, exponent = e }+      | c == 0    -> (0, 0)+      | otherwise -> (m, n)+      where m = signFunc s (fromIntegral $ c * 10^d)+            n = fromIntegral e - d+            d = floatDigits x - numDigits c+    Inf { sign = s                                } -> (special s 0, maxBound)+    NaN { sign = s, signaling = sig, payload = p  } -> (special s p, n)+      where n = minBound + fromEnum sig+     where special :: Sign -> Coefficient -> Integer           special s v = signFunc s (pp + fromIntegral v)           pp = 10 ^ floatDigits x :: Integer@@ -378,16 +517,16 @@                                    }           special             | n == maxBound     = Inf  { sign = signMatch m }-            | n == minBound     = QNaN { sign = signMatch m, payload = p }-            | otherwise         = SNaN { sign = signMatch m, payload = p }+            | n == minBound     = qNaN { sign = signMatch m, payload = p }+            | otherwise         = sNaN { sign = signMatch m, payload = p }             where p = fromInteger (am - pp)+           am = abs m              :: Integer           pp = 10 ^ floatDigits x :: Integer    isNaN x = case x of-    QNaN{} -> True-    SNaN{} -> True-    _      -> False+    NaN{} -> True+    _     -> False    isInfinite x = case x of     Inf{} -> True@@ -401,17 +540,6 @@    isIEEE _ = True -{- $doctest-RealFloat-prop> uncurry encodeFloat (decodeFloat x) == (x :: BasicDecimal)-prop> isFinite x ==> significand x * fromInteger (floatRadix x) ^^ Prelude.exponent x == (x :: BasicDecimal)-prop> Prelude.exponent (0 :: BasicDecimal) == 0-prop> isFinite x && x /= 0 ==> Prelude.exponent (x :: BasicDecimal) == snd (decodeFloat x) + floatDigits x--prop> isNegativeZero (read "-0" :: BasicDecimal) == True-prop> isNegativeZero (read "+0" :: BasicDecimal) == False-prop> x /= 0 ==> isNegativeZero (x :: BasicDecimal) == False--}- -- | The 'Bits' instance makes use of the logical operations from the -- /General Decimal Arithmetic Specification/ using a /digit-wise/ -- representation of bits where the /sign/ is non-negative, the /exponent/ is@@ -454,6 +582,13 @@ instance FinitePrecision p => FiniteBits (Decimal p r) where   finiteBitSize x = let Just p = precision x in p +instance NFData (Decimal p r) where+  rnf Num { sign = s, coefficient = c, exponent = e } =+    rnf s `seq` rnf c `seq` rnf e+  rnf Inf { sign = s } = rnf s+  rnf NaN { sign = s, signaling = sig, payload = p } =+    rnf s `seq` rnf sig `seq` rnf p+ -- | A 'Decimal' representing the value zero zero :: Decimal p r zero = Num { sign        = Pos@@ -465,6 +600,10 @@ oneHalf :: Decimal p r oneHalf = zero { coefficient = 5, exponent = -1 } +-- | A 'Decimal' representing the value ¼+oneQuarter :: Decimal p r+oneQuarter = zero { coefficient = 25, exponent = -2 }+ -- | A 'Decimal' representing the value one one :: Decimal p r one = zero { coefficient = 1 }@@ -473,6 +612,10 @@ two :: Decimal p r two = zero { coefficient = 2 } +-- | A 'Decimal' representing the value three+three :: Decimal p r+three = zero { coefficient = 3 }+ -- | A 'Decimal' representing the value ten ten :: Decimal p r ten = zero { coefficient = 10 }@@ -487,11 +630,11 @@  -- | A 'Decimal' representing undefined results qNaN :: Decimal p r-qNaN = QNaN { sign = Pos, payload = 0 }+qNaN = NaN { sign = Pos, signaling = False, payload = 0 }  -- | A signaling 'Decimal' representing undefined results sNaN :: Decimal p r-sNaN = SNaN { sign = Pos, payload = 0 }+sNaN = qNaN { signaling = True }  -- | Negate the given 'Decimal' by directly flipping its sign. flipSign :: Decimal p r -> Decimal p r@@ -559,35 +702,31 @@   | isFinite n && not (isZero n) = maybe False (adjustedExponent n <) (eMin n)   | otherwise                    = False --- | If the argument is 'False', return a 'Decimal' value zero; if 'True',--- return the value one. This is basically an optimized @toEnum . fromEnum@ to--- support an all-decimal usage of the operations from--- "Numeric.Decimal.Operation" that return a 'Bool'.+-- | Return @0@ or @1@ if the argument is 'False' or 'True', respectively.+-- This is basically an optimized @'toEnum' . 'fromEnum'@ and allows an+-- all-decimal usage of the operations from "Numeric.Decimal.Operation" that+-- return a 'Bool'. fromBool :: Bool -> Decimal p r fromBool False = zero fromBool True  = one  -- | Return 'False' if the argument is zero or NaN, and 'True' otherwise. toBool :: Decimal p r -> Bool-toBool Num { coefficient = c }-  | c == 0    = False-  | otherwise = True-toBool Inf{}  = True-toBool _      = False+toBool Num { coefficient = c } = c /= 0+toBool NaN{}                   = False+toBool _                       = True +-- | Return @-1@, @0@, or @1@ if the argument is 'LT', 'EQ', or 'GT',+-- respectively. This allows an all-decimal usage of the operations from+-- "Numeric.Decimal.Operation" that return an 'Ordering'.+fromOrdering :: Ordering -> Decimal p r+fromOrdering LT = negativeOne+fromOrdering EQ = zero+fromOrdering GT = one+ -- | Upper limit on the absolute value of the exponent eLimit :: Precision p => p -> Maybe Exponent eLimit = eMax -- ?---- | Minimum value of the adjusted exponent-eMin :: Precision p => p -> Maybe Exponent-eMin n = (1 -) <$> eMax n---- | Maximum value of the adjusted exponent-eMax :: Precision p => p -> Maybe Exponent-eMax n = subtract 1 . (10 ^) . numDigits <$> base-  where mlength = precision n                    :: Maybe Int-        base = (10 *) . fromIntegral <$> mlength :: Maybe Coefficient  -- | Minimum value of the exponent for subnormal results eTiny :: Precision p => p -> Maybe Exponent

@@ -4,37 +4,31 @@  module Numeric.Decimal.Number        ( Sign(..)-       , Decimal(..)        , Coefficient+       , Exponent+       , Decimal(..)        , numDigits        ) where  import Numeric.Natural (Natural) -import Numeric.Decimal.Precision- data Sign = Pos | Neg-instance Eq Sign  type Coefficient = Natural-type Exponent    = Int+type Exponent    = Integer type Payload     = Coefficient  type role Decimal phantom phantom data Decimal p r-  = Num  { sign        :: Sign-         , coefficient :: Coefficient-         , exponent    :: Exponent-         }-  | Inf  { sign        :: Sign-         }-  | QNaN { sign        :: Sign-         , payload     :: Payload-         }-  | SNaN { sign        :: Sign-         , payload     :: Payload-         }--instance Precision p => Precision (Decimal p r)+  = Num { sign        :: Sign+        , coefficient :: Coefficient+        , exponent    :: Exponent+        }+  | Inf { sign        :: Sign+        }+  | NaN { sign        :: Sign+        , signaling   :: Bool+        , payload     :: Payload+        }  numDigits :: Coefficient -> Int

@@ -16,1950 +16,1383 @@        ( -- * Arithmetic operations          -- $arithmetic-operations -         abs-       , add-       , subtract-       , compare-       , compareSignal-       , divide-         -- divideInteger-       , exp-       , fusedMultiplyAdd-       , ln-       , log10-       , max-       , maxMagnitude-       , min-       , minMagnitude-       , minus-       , plus-       , multiply-         -- nextMinus-         -- nextPlus-         -- nextToward-       , power-       , quantize-       , reduce-         -- remainder-         -- remainderNear-         -- roundToIntegralExact-         -- roundToIntegralValue-       , squareRoot--         -- * Miscellaneous operations-         -- $miscellaneous-operations--       , and-       , canonical-       , class_, Class(..), Sign(..), NumberClass(..), NaNClass(..)-         -- compareTotal-         -- compareTotalMagnitude-       , copy-       , copyAbs-       , copyNegate-       , copySign-       , invert-       , isCanonical-       , isFinite-       , isInfinite-       , isNaN-       , isNormal-       , isQNaN-       , isSigned-       , isSNaN-       , isSubnormal-       , isZero-       , logb-       , or-       , radix-       , rotate-       , sameQuantum-         -- scaleb-       , shift-       , xor-       ) where--import Prelude hiding (abs, and, compare, exp, exponent, isInfinite, isNaN,-                       max, min, or, subtract)-import qualified Prelude--import Control.Monad (join)-import Data.Bits (complement, setBit, testBit, zeroBits, (.&.), (.|.))-import Data.Coerce (coerce)-import Data.List (find)-import Data.Maybe (fromMaybe)--import qualified Data.Bits as Bits--import Numeric.Decimal.Arithmetic-import Numeric.Decimal.Number hiding (isFinite, isNormal, isSubnormal, isZero)-import Numeric.Decimal.Precision-import Numeric.Decimal.Rounding--import qualified Numeric.Decimal.Number as Number--{- $setup->>> :load Harness--}--finitePrecision :: FinitePrecision p => Decimal p r -> Int-finitePrecision n = let Just p = precision n in p--roundingAlg :: Rounding r => Arith p r a -> RoundingAlgorithm-roundingAlg = rounding . arithRounding-  where arithRounding :: Arith p r a -> r-        arithRounding = undefined--result :: (Precision p, Rounding r) => Decimal p r -> Arith p r (Decimal p r)-result = roundDecimal  -- ...---  | maybe False (numDigits c >) (precision r) = undefined--invalidOperation :: Decimal a b -> Arith p r (Decimal p r)-invalidOperation n = raiseSignal InvalidOperation qNaN--toQNaN :: Decimal a b -> Decimal p r-toQNaN SNaN { sign = s, payload = p } = QNaN { sign = s, payload = p }-toQNaN n@QNaN{}                       = coerce n-toQNaN n                              = qNaN { sign = sign n }--toQNaN2 :: Decimal a b -> Decimal c d -> Decimal p r-toQNaN2 nan@SNaN{} _ = toQNaN nan-toQNaN2 _ nan@SNaN{} = toQNaN nan-toQNaN2 nan@QNaN{} _ = coerce nan-toQNaN2 _ nan@QNaN{} = coerce nan-toQNaN2 n _          = toQNaN n--quietToSignal :: Decimal p r -> Decimal p r-quietToSignal QNaN { sign = s, payload = p } = SNaN { sign = s, payload = p }-quietToSignal x = x---- $arithmetic-operations------ This section describes the arithmetic operations on, and some other--- functions of, numbers, including subnormal numbers, negative zeros, and--- special values (see also IEEE 754 §5 and §6).--{- $doctest-special-values->>> op2 Op.add "Infinity" "1"-Infinity-->>> op2 Op.add "NaN" "1"-NaN-->>> op2 Op.add "NaN" "Infinity"-NaN-->>> op2 Op.subtract "1" "Infinity"--Infinity-->>> op2 Op.multiply "-1" "Infinity"--Infinity-->>> op2 Op.subtract "-0" "0"--0-->>> op2 Op.multiply "-1" "0"--0-->>> op2 Op.divide "1" "0"-Infinity-->>> op2 Op.divide "1" "-0"--Infinity-->>> op2 Op.divide "-1" "0"--Infinity--}---- | 'add' takes two operands. If either operand is a /special value/ then the--- general rules apply.------ Otherwise, the operands are added.------ The result is then rounded to /precision/ digits if necessary, counting--- from the most significant digit of the result.-add :: (Precision p, Rounding r)-    => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)-add Num { sign = xs, coefficient = xc, exponent = xe }-    Num { sign = ys, coefficient = yc, exponent = ye } = sum--  where sum = result Num { sign = rs, coefficient = rc, exponent = re }-        rs | rc /= 0                       = if xac > yac then xs else ys-           | xs == Neg && ys == Neg        = Neg-           | xs /= ys &&-             roundingAlg sum == RoundFloor = Neg-           | otherwise                     = Pos-        rc | xs == ys  = xac + yac-           | xac > yac = xac - yac-           | otherwise = yac - xac-        re = Prelude.min xe ye-        (xac, yac) | xe == ye  = (xc, yc)-                   | xe >  ye  = (xc * 10^n, yc)-                   | otherwise = (xc, yc * 10^n)-          where n = Prelude.abs (xe - ye)--add inf@Inf { sign = xs } Inf { sign = ys }-  | xs == ys  = return (coerce inf)-  | otherwise = invalidOperation inf-add inf@Inf{} Num{} = return (coerce inf)-add Num{} inf@Inf{} = return (coerce inf)-add x y             = return (toQNaN2 x y)--{- $doctest-add->>> op2 Op.add "12" "7.00"-19.00-->>> op2 Op.add "1E+2" "1E+4"-1.01E+4--}---- | 'subtract' takes two operands. If either operand is a /special value/--- then the general rules apply.------ Otherwise, the operands are added after inverting the /sign/ used for the--- second operand.------ The result is then rounded to /precision/ digits if necessary, counting--- from the most significant digit of the result.-subtract :: (Precision p, Rounding r)-         => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)-subtract x = add x . flipSign--{- $doctest-subtract->>> op2 Op.subtract "1.3" "1.07"-0.23-->>> op2 Op.subtract "1.3" "1.30"-0.00-->>> op2 Op.subtract "1.3" "2.07"--0.77--}---- | 'minus' takes one operand, and corresponds to the prefix minus operator--- in programming languages.------ Note that the result of this operation is affected by context and may set--- /flags/. The 'copyNegate' operation may be used instead of 'minus' if this--- is not desired.-minus :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r)-minus x = zero { exponent = exponent x } `subtract` x--{- $doctest-minus->>> op1 Op.minus "1.3"--1.3-->>> op1 Op.minus "-1.3"-1.3--}---- | 'plus' takes one operand, and corresponds to the prefix plus operator in--- programming languages.------ Note that the result of this operation is affected by context and may set--- /flags/.-plus :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r)-plus x = zero { exponent = exponent x } `add` x--{- $doctest-plus->>> op1 Op.plus "1.3"-1.3-->>> op1 Op.plus "-1.3"--1.3--}---- | 'multiply' takes two operands. If either operand is a /special value/--- then the general rules apply. Otherwise, the operands are multiplied--- together (“long multiplication”), resulting in a number which may be as--- long as the sum of the lengths of the two operands.------ The result is then rounded to /precision/ digits if necessary, counting--- from the most significant digit of the result.-multiply :: (Precision p, Rounding r)-         => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)-multiply Num { sign = xs, coefficient = xc, exponent = xe }-         Num { sign = ys, coefficient = yc, exponent = ye } = result rn--  where rn = Num { sign = rs, coefficient = rc, exponent = re }-        rs = xorSigns xs ys-        rc = xc * yc-        re = xe + ye--multiply Inf { sign = xs } Inf { sign = ys } =-  return Inf { sign = xorSigns xs ys }-multiply Inf { sign = xs } Num { sign = ys, coefficient = yc }-  | yc == 0   = invalidOperation qNaN-  | otherwise = return Inf { sign = xorSigns xs ys }-multiply Num { sign = xs, coefficient = xc } Inf { sign = ys }-  | xc == 0   = invalidOperation qNaN-  | otherwise = return Inf { sign = xorSigns xs ys }-multiply nan@SNaN{} _ = invalidOperation nan-multiply _ nan@SNaN{} = invalidOperation nan-multiply x y = return (toQNaN2 x y)--{- $doctest-multiply->>> op2 Op.multiply "1.20" "3"-3.60-->>> op2 Op.multiply "7" "3"-21-->>> op2 Op.multiply "0.9" "0.8"-0.72-->>> op2 Op.multiply "0.9" "-0"--0.0-->>> op2 Op.multiply "654321" "654321"-4.28135971E+11--}---- | 'exp' takes one operand. If the operand is a NaN then the general rules--- for special values apply.------ Otherwise, the result is /e/ raised to the power of the operand, with the--- following cases:------ * If the operand is −Infinity, the result is 0 and exact.------ * If the operand is a zero, the result is 1 and exact.------ * If the operand is +Infinity, the result is +Infinity and exact.------ * Otherwise the result is inexact and will be rounded using the--- /round-half-even/ algorithm. The coefficient will have exactly /precision/--- digits (unless the result is subnormal). These inexact results should be--- correctly rounded, but may be up to 1 ulp (unit in last place) in error.-exp :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p RoundHalfEven)-exp x@Num { sign = s, coefficient = c }-  | c == 0    = return one-  | s == Neg  = subArith (maclaurin x { sign = Pos } >>= reciprocal) >>=-                subRounded >>= result-  | otherwise = subArith (maclaurin x) >>= subRounded >>= result--  where multiplyExact :: Decimal a b -> Decimal c d-                      -> Arith PInfinite RoundHalfEven-                         (Decimal PInfinite RoundHalfEven)-        multiplyExact = multiply--        maclaurin :: FinitePrecision p => Decimal a b-                  -> Arith p RoundHalfEven (Decimal p RoundHalfEven)-        maclaurin x-          | adjustedExponent x >= 0 = subArith (subMaclaurin x) >>= subRounded-          | otherwise = sum one one one one-          where sum :: FinitePrecision p-                    => Decimal p RoundHalfEven-                    -> Decimal PInfinite RoundHalfEven-                    -> Decimal PInfinite RoundHalfEven-                    -> Decimal PInfinite RoundHalfEven-                    -> Arith p RoundHalfEven (Decimal p RoundHalfEven)-                sum s num den n = do-                  num' <- subArith (multiplyExact num x)-                  den' <- subArith (multiplyExact den n)-                  s' <- add s =<< divide num' den'-                  if s' == s then return s'-                    else sum s' num' den' =<< subArith (add n one)--        subMaclaurin :: FinitePrecision p => Decimal a b-                     -> Arith p RoundHalfEven (Decimal p RoundHalfEven)-        subMaclaurin x = subArith (multiplyExact x oneHalf) >>= maclaurin >>=-          \r -> multiply r r--        subRounded :: Precision p-                   => Decimal (PPlus1 (PPlus1 p)) a-                   -> Arith p r (Decimal p RoundHalfEven)-        subRounded = subArith . roundDecimal--        result :: Decimal p a -> Arith p r (Decimal p a)-        result r = coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r')-          where r' = coerce r--exp n@Inf { sign = s }-  | s == Pos  = return (coerce n)-  | otherwise = return zero-exp n@QNaN{}  = return (coerce n)-exp n@SNaN{}  = coerce <$> invalidOperation n--{- $doctest-exp->>> op1 Op.exp "-Infinity"-0-->>> op1 Op.exp "-1"-0.367879441-->>> op1 Op.exp "0"-1-->>> op1 Op.exp "1"-2.71828183-->>> op1 Op.exp "0.693147181"-2.00000000-->>> op1 Op.exp "+Infinity"-Infinity--}---- | 'fusedMultiplyAdd' takes three operands; the first two are multiplied--- together, using 'multiply', with sufficient precision and exponent range--- that the result is exact and unrounded. No /flags/ are set by the--- multiplication unless one of the first two operands is a signaling NaN or--- one is a zero and the other is an infinity.------ Unless the multiplication failed, the third operand is then added to the--- result of that multiplication, using 'add', under the current context.------ In other words, @fusedMultiplyAdd x y z@ delivers a result which is @(x ×--- y) + z@ with only the one, final, rounding.-fusedMultiplyAdd :: (Precision p, Rounding r)-                 => Decimal a b -> Decimal c d -> Decimal e f-                 -> Arith p r (Decimal p r)-fusedMultiplyAdd x y z =-  either raise (return . coerce) (exactMult x y) >>= add z--  where exactMult :: Rounding r => Decimal a b -> Decimal c d-                  -> Either (Exception PInfinite r) (Decimal PInfinite r)-        exactMult x y = evalArith (multiply x y) newContext--        raise :: Exception a r -> Arith p r (Decimal p r)-        raise e = raiseSignal (exceptionSignal e) (coerce $ exceptionResult e)--{- $doctest-fusedMultiplyAdd->>> op3 Op.fusedMultiplyAdd "3" "5" "7"-22-->>> op3 Op.fusedMultiplyAdd "3" "-5" "7"--8-->>> op3 Op.fusedMultiplyAdd "888565290" "1557.96930" "-86087.7578"-1.38435736E+12--}---- | 'ln' takes one operand. If the operand is a NaN then the general rules--- for special values apply.------ Otherwise, the operand must be a zero or positive, and the result is the--- natural (base /e/) logarithm of the operand, with the following cases:------ * If the operand is a zero, the result is −Infinity and exact.------ * If the operand is +Infinity, the result is +Infinity and exact.------ * If the operand equals one, the result is 0 and exact.------ * Otherwise the result is inexact and will be rounded using the--- /round-half-even/ algorithm. The coefficient will have exactly /precision/--- digits (unless the result is subnormal). These inexact results should be--- correctly rounded, but may be up to 1 ulp (unit in last place) in error.-ln :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p RoundHalfEven)-ln x@Num { sign = s, coefficient = c, exponent = e }-  | c == 0   = return infinity { sign = Neg }-  | s == Pos = if e <= 0 && c == 10^(-e) then return zero-               else subArith (subLn x) >>= subRounded >>= result--  where subLn :: FinitePrecision p => Decimal a b-              -> Arith p RoundHalfEven (Decimal p RoundHalfEven)-        subLn x = do-          let fe = fromIntegral (-(numDigits c - 1)) :: Exponent-              r  = fromIntegral (e - fe) :: Decimal PInfinite RoundHalfEven-          lnf <- taylorLn x { exponent = fe }-          add lnf =<< multiply r =<< ln10--        subRounded :: Precision p => Decimal (PPlus1 (PPlus1 p)) a-                   -> Arith p r (Decimal p RoundHalfEven)-        subRounded = subArith . roundDecimal--        result :: Decimal p a -> Arith p r (Decimal p a)-        result r = coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r')-          where r' = coerce r--ln n@Inf { sign = Pos } = return (coerce n)-ln n@QNaN{} = return (coerce n)-ln n = coerce <$> invalidOperation n--taylorLn :: FinitePrecision p => Decimal a b-         -> Arith p RoundHalfEven (Decimal p RoundHalfEven)-taylorLn x = do-  num <- x `subtract` one-  den <- x `add`      one-  multiply two =<< sum =<< num `divide` den--    where sum :: FinitePrecision p => Decimal p RoundHalfEven-              -> Arith p RoundHalfEven (Decimal p RoundHalfEven)-          sum b = multiply b b >>= \b2 -> sum' b b b2 one--            where sum' :: FinitePrecision p-                       => Decimal p RoundHalfEven-                       -> Decimal p RoundHalfEven-                       -> Decimal p RoundHalfEven-                       -> Decimal PInfinite RoundHalfEven-                       -> Arith p RoundHalfEven (Decimal p RoundHalfEven)-                  sum' s m b n = do-                    m' <- multiply m b-                    n' <- subArith (add n two)-                    s' <- add s =<< divide m' n'-                    if s' == s then return s' else sum' s' m' b n'--ln10 :: FinitePrecision p => Arith p r (Decimal p RoundHalfEven)-ln10 = getPrecision >>= \(Just p) ->-  if p <= 50 then return fastLn10 else slowLn10--  where fastLn10 :: FinitePrecision p => Decimal p RoundHalfEven-        fastLn10 = 2.3025850929940456840179914546843642076011014886288--        slowLn10 :: FinitePrecision p => Arith p r (Decimal p RoundHalfEven)-        slowLn10 = subArith (taylorLn ten) >>= subRound--          where subRound :: Precision p => Decimal (PPlus1 (PPlus1 p)) a-                         -> Arith p r (Decimal p RoundHalfEven)-                subRound = subArith . roundDecimal--{- $doctest-ln->>> op1 Op.ln "0"--Infinity-->>> op1 Op.ln "1.000"-0-->>> op1 Op.ln "2.71828183"-1.00000000-->>> op1 Op.ln "10"-2.30258509-->>> op1 Op.ln "+Infinity"-Infinity--}---- | 'log10' takes one operand. If the operand is a NaN then the general rules--- for special values apply.------ Otherwise, the operand must be a zero or positive, and the result is the--- base 10 logarithm of the operand, with the following cases:------ * If the operand is a zero, the result is −Infinity and exact.------ * If the operand is +Infinity, the result is +Infinity and exact.------ * If the operand equals an integral power of ten (including 10^0 and--- negative powers) and there is sufficient /precision/ to hold the integral--- part of the result, the result is an integer (with an exponent of 0) and--- exact.------ * Otherwise the result is inexact and will be rounded using the--- /round-half-even/ algorithm. The coefficient will have exactly /precision/--- digits (unless the result is subnormal). These inexact results should be--- correctly rounded, but may be up to 1 ulp (unit in last place) in error.-log10 :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p RoundHalfEven)-log10 x@Num { sign = s, coefficient = c, exponent = e }-  | c == 0   = return infinity { sign = Neg }-  | s == Pos = getPrecision >>= \prec -> case powerOfTen c e of-      Just p | maybe True (numDigits pc <=) prec -> return (fromInteger p)-        where pc = fromInteger (Prelude.abs p) :: Coefficient-      _ -> subArith (join $ divide <$> ln x <*> ln10) >>= result--  where powerOfTen :: Coefficient -> Exponent -> Maybe Integer-        powerOfTen c e-          | c == 10^d = Just (fromIntegral e + fromIntegral d)-          | otherwise = Nothing-          where d = numDigits c - 1 :: Int--        result :: Decimal p a -> Arith p r (Decimal p a)-        result r = coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r')-          where r' = coerce r--log10 n@Inf { sign = Pos } = return (coerce n)-log10 n@QNaN{} = return (coerce n)-log10 n = coerce <$> invalidOperation n--{- $doctest-log10->>> op1 Op.log10 "0"--Infinity-->>> op1 Op.log10 "0.001"--3-->>> op1 Op.log10 "1.000"-0-->>> op1 Op.log10 "2"-0.301029996-->>> op1 Op.log10 "10"-1-->>> op1 Op.log10 "70"-1.84509804-->>> op1 Op.log10 "+Infinity"-Infinity--}---- | 'divide' takes two operands. If either operand is a /special value/ then--- the general rules apply.------ Otherwise, if the divisor is zero then either the Division undefined--- condition is raised (if the dividend is zero) and the result is NaN, or the--- Division by zero condition is raised and the result is an Infinity with a--- sign which is the exclusive or of the signs of the operands.------ Otherwise, a “long division” is effected.------ The result is then rounded to /precision/ digits, if necessary, according--- to the /rounding/ algorithm and taking into account the remainder from the--- division.-divide :: (FinitePrecision p, Rounding r)-       => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)-divide dividend@Num{ sign = xs } Num { coefficient = 0, sign = ys }-  | Number.isZero dividend = invalidOperation qNaN-  | otherwise              = raiseSignal DivisionByZero-                             infinity { sign = xorSigns xs ys }-divide Num { sign = xs, coefficient = xc, exponent = xe }-       Num { sign = ys, coefficient = yc, exponent = ye } = quotient--  where quotient = result =<< answer-        rn = Num { sign = rs, coefficient = rc, exponent = re }-        rs = xorSigns xs ys-        (rc, rem, dv, adjust) = longDivision xc yc (finitePrecision rn)-        re = xe - (ye + adjust)-        answer-          | rem == 0  = return rn-          | otherwise = roundDecimal $ case (rem * 2) `Prelude.compare` dv of-              LT -> rn { coefficient = rc * 10 + 1, exponent = re - 1 }-              EQ -> rn { coefficient = rc * 10 + 5, exponent = re - 1 }-              GT -> rn { coefficient = rc * 10 + 9, exponent = re - 1 }--divide Inf{} Inf{} = invalidOperation qNaN-divide Inf { sign = xs } Num { sign = ys } =-  return Inf { sign = xorSigns xs ys }-divide Num { sign = xs } Inf { sign = ys } =-  return zero { sign = xorSigns xs ys }-divide x y = return (toQNaN2 x y)--{- $doctest-divide->>> op2 Op.divide "1" "3"-0.333333333-->>> op2 Op.divide "2" "3"-0.666666667-->>> op2 Op.divide "5" "2"-2.5-->>> op2 Op.divide "1" "10"-0.1-->>> op2 Op.divide "12" "12"-1-->>> op2 Op.divide "8.00" "2"-4.00-->>> op2 Op.divide "2.400" "2.0"-1.20-->>> op2 Op.divide "1000" "100"-10-->>> op2 Op.divide "1000" "1"-1000-->>> op2 Op.divide "2.40E+6" "2"-1.20E+6--}--type Dividend  = Coefficient-type Divisor   = Coefficient-type Quotient  = Coefficient-type Remainder = Dividend--longDivision :: Dividend -> Divisor -> Int-             -> (Quotient, Remainder, Divisor, Exponent)-longDivision 0  dv _ = (0, 0, dv, 0)-longDivision dd dv p = step1 dd dv 0--  where step1 :: Dividend -> Divisor -> Exponent-              -> (Quotient, Remainder, Divisor, Exponent)-        step1 dd dv adjust-          | dd <       dv = step1 (dd * 10)  dv       (adjust + 1)-          | dd >= 10 * dv = step1  dd       (dv * 10) (adjust - 1)-          | otherwise     = step2  dd        dv        adjust--        step2 :: Dividend -> Divisor -> Exponent-              -> (Quotient, Remainder, Divisor, Exponent)-        step2 = step3 0--        step3 :: Quotient -> Dividend -> Divisor -> Exponent-              -> (Quotient, Remainder, Divisor, Exponent)-        step3 r dd dv adjust-          | dv <= dd                 = step3 (r +  1) (dd - dv) dv  adjust-          | (dd == 0 && adjust >= 0) ||-            numDigits r == p         = step4  r        dd       dv  adjust-          | otherwise                = step3 (r * 10) (dd * 10) dv (adjust + 1)--        step4 :: Quotient -> Remainder -> Divisor -> Exponent-              -> (Quotient, Remainder, Divisor, Exponent)-        step4 = (,,,)--reciprocal :: (FinitePrecision p, Rounding r)-           => Decimal a b -> Arith p r (Decimal p r)-reciprocal = divide one---- | 'abs' takes one operand. If the operand is negative, the result is the--- same as using the 'minus' operation on the operand. Otherwise, the result--- is the same as using the 'plus' operation on the operand.------ Note that the result of this operation is affected by context and may set--- /flags/. The 'copyAbs' operation may be used if this is not desired.-abs :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r)-abs x-  | isNegative x = minus x-  | otherwise    = plus  x--{- $doctest-abs->>> op1 Op.abs "2.1"-2.1-->>> op1 Op.abs "-100"-100-->>> op1 Op.abs "101.5"-101.5-->>> op1 Op.abs "-101.5"-101.5--}---- | 'compare' takes two operands and compares their values numerically. If--- either operand is a /special value/ then the general rules apply. No flags--- are set unless an operand is a signaling NaN.------ Otherwise, the operands are compared, returning @−1@ if the first is less--- than the second, @0@ if they are equal, or @1@ if the first is greater than--- the second.-compare :: (Precision p, Rounding r)-        => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)-compare x@Num{} y@Num{} = nzp <$> (xn `subtract` yn)--  where (xn, yn) | sign x /= sign y = (nzp x, nzp y)-                 | otherwise        = (x, y)--        nzp :: Decimal p r -> Decimal p r-        nzp Num { sign = s, coefficient = c }-          | c == 0    = zero-          | s == Pos  = one-          | otherwise = negativeOne-        nzp Inf { sign = s }-          | s == Pos  = one-          | otherwise = negativeOne-        nzp n = toQNaN n--compare Inf { sign = xs } Inf { sign = ys }-  | xs == ys  = return zero-  | xs == Neg = return negativeOne-  | otherwise = return one-compare Inf { sign = xs } Num { }-  | xs == Neg = return negativeOne-  | otherwise = return one-compare Num { } Inf { sign = ys }-  | ys == Pos = return negativeOne-  | otherwise = return one-compare nan@SNaN{} _ = invalidOperation nan-compare _ nan@SNaN{} = invalidOperation nan-compare x y          = return (toQNaN2 x y)--{- $doctest-compare->>> op2 Op.compare "2.1" "3"--1-->>> op2 Op.compare "2.1" "2.1"-0-->>> op2 Op.compare "2.1" "2.10"-0-->>> op2 Op.compare "3" "2.1"-1-->>> op2 Op.compare "2.1" "-3"-1-->>> op2 Op.compare "-3" "2.1"--1--}---- | 'compareSignal' takes two operands and compares their values--- numerically. This operation is identical to 'compare', except that if--- neither operand is a signaling NaN then any quiet NaN operand is treated as--- though it were a signaling NaN. (That is, all NaNs signal, with signaling--- NaNs taking precedence over quiet NaNs.)-compareSignal :: (Precision p, Rounding r)-              => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)-compareSignal x@SNaN{} y        =               x `compare`               y-compareSignal x        y@SNaN{} =               x `compare`               y-compareSignal x        y        = quietToSignal x `compare` quietToSignal y---- | 'max' takes two operands, compares their values numerically, and returns--- the maximum. If either operand is a NaN then the general rules apply,--- unless one is a quiet NaN and the other is numeric, in which case the--- numeric operand is returned.-max :: (Precision p, Rounding r)-    => Decimal a b -> Decimal a b -> Arith p r (Decimal a b)-max x y = snd <$> minMax id x y--{- $doctest-max->>> op2 Op.max "3" "2"-3-->>> op2 Op.max "-10" "3"-3-->>> op2 Op.max "1.0" "1"-1-->>> op2 Op.max "7" "NaN"-7--}---- | 'maxMagnitude' takes two operands and compares their values numerically--- with their /sign/ ignored and assumed to be 0.------ If, without signs, the first operand is the larger then the original first--- operand is returned (that is, with the original sign). If, without signs,--- the second operand is the larger then the original second operand is--- returned. Otherwise the result is the same as from the 'max' operation.-maxMagnitude :: (Precision p, Rounding r)-             => Decimal a b -> Decimal a b -> Arith p r (Decimal a b)-maxMagnitude x y = snd <$> minMax withoutSign x y---- | 'min' takes two operands, compares their values numerically, and returns--- the minimum. If either operand is a NaN then the general rules apply,--- unless one is a quiet NaN and the other is numeric, in which case the--- numeric operand is returned.-min :: (Precision p, Rounding r)-    => Decimal a b -> Decimal a b -> Arith p r (Decimal a b)-min x y = fst <$> minMax id x y--{- $doctest-min->>> op2 Op.min "3" "2"-2-->>> op2 Op.min "-10" "3"--10-->>> op2 Op.min "1.0" "1"-1.0-->>> op2 Op.min "7" "NaN"-7--}---- | 'minMagnitude' takes two operands and compares their values numerically--- with their /sign/ ignored and assumed to be 0.------ If, without signs, the first operand is the smaller then the original first--- operand is returned (that is, with the original sign). If, without signs,--- the second operand is the smaller then the original second operand is--- returned. Otherwise the result is the same as from the 'min' operation.-minMagnitude :: (Precision p, Rounding r)-             => Decimal a b -> Decimal a b -> Arith p r (Decimal a b)-minMagnitude x y = fst <$> minMax withoutSign x y---- | Ordering function for 'min', 'minMagnitude', 'max', and 'maxMagnitude':--- returns the original arguments as (smaller, larger) when the given function--- is applied to them.-minMax :: (Precision p, Rounding r)-       => (Decimal a b -> Decimal a b) -> Decimal a b -> Decimal a b-       -> Arith p r (Decimal a b, Decimal a b)-minMax _ x@Num{}  QNaN{} = return (x, x)-minMax _ x@Inf{}  QNaN{} = return (x, x)-minMax _  QNaN{} y@Num{} = return (y, y)-minMax _  QNaN{} y@Inf{} = return (y, y)--minMax f x y = do-  c <- f x `compare` f y-  return $ case c of-    Num { coefficient = 0 } -> case (sign x, sign y) of-      (Neg, Pos) -> (x, y)-      (Pos, Neg) -> (y, x)-      (Pos, Pos) -> case (x, y) of-        (Num { exponent = xe }, Num { exponent = ye }) | xe > ye -> (y, x)-        _ -> (x, y)-      (Neg, Neg) -> case (x, y) of-        (Num { exponent = xe }, Num { exponent = ye }) | xe < ye -> (y, x)-        _ -> (x, y)-    Num { sign = Pos } -> (y, x)-    Num { sign = Neg } -> (x, y)-    nan -> let nan' = coerce nan in (nan', nan')--withoutSign :: Decimal p r -> Decimal p r-withoutSign n = n { sign = Pos }---- | 'power' takes two operands, and raises a number (the left-hand operand)--- to a power (the right-hand operand). If either operand is a /special value/--- then the general rules apply, except in certain cases.-power :: (FinitePrecision p, Rounding r)-      => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)-power x@Num { coefficient = 0 } y@Num{}-  | Number.isZero y     = invalidOperation qNaN-  | Number.isNegative y = return infinity { sign = powerSign x y }-  | otherwise           = return zero     { sign = powerSign x y }-power x@Num{} y@Num{} = case integralValue y of-  Just i  | i < 0               -> reciprocal x >>= \rx -> integralPower rx (-i)-          | otherwise           ->                         integralPower  x   i-  Nothing | Number.isPositive x -> ln x >>= multiply y >>= fmap coerce . exp-          | otherwise           -> invalidOperation qNaN-power x@Num{} y@Inf{}-  | Number.isPositive x = return $ case sign y of-      Pos -> infinity-      Neg -> zero-  | otherwise           = invalidOperation qNaN-power x@Inf{} y@Num{}-  | Number.isZero y     = return one-  | Number.isPositive y = return infinity { sign = powerSign x y }-  | otherwise           = return zero     { sign = powerSign x y }-power Inf{} Inf { sign = s }-  | s == Pos            = return infinity-  | otherwise           = return zero-power x@SNaN{} _        = invalidOperation x-power _        y@SNaN{} = invalidOperation y-power x@QNaN{} _        = return (coerce x)-power _        y@QNaN{} = return (coerce y)--powerSign :: Decimal a b -> Decimal c d -> Sign-powerSign x y-  | Number.isNegative x && fromMaybe False (odd <$> integralValue y) = Neg-  | otherwise                                                        = Pos--integralPower :: (Precision p, Rounding r)-              => Decimal a b -> Integer -> Arith p r (Decimal p r)-integralPower b e = integralPower' (return b) e one-  where integralPower' :: (Precision p, Rounding r)-                       => Arith p r (Decimal a b) -> Integer -> Decimal p r-                       -> Arith p r (Decimal p r)-        integralPower' _  0 r = return r-        integralPower' mb e r-          | odd e     = mb >>= \b -> multiply r b >>=-                        integralPower'              (multiply b b) e'-          | otherwise = integralPower' (mb >>= \b -> multiply b b) e' r-          where e' = e `div` 2--{- $doctest-power->>> op2 Op.power "2" "3"-8-->>> op2 Op.power "-2" "3"--8-->>> op2 Op.power "2" "-3"-0.125-->>> op2 Op.power "1.7" "8"-69.7575744-->>> op2 Op.power "10" "0.301029996"-2.00000000-->>> op2 Op.power "Infinity" "-1"-0-->>> op2 Op.power "Infinity" "0"-1-->>> op2 Op.power "Infinity" "1"-Infinity-->>> op2 Op.power "-Infinity" "-1"--0-->>> op2 Op.power "-Infinity" "0"-1-->>> op2 Op.power "-Infinity" "1"--Infinity-->>> op2 Op.power "-Infinity" "2"-Infinity-->>> op2 Op.power "0" "0"-NaN--}---- | 'quantize' takes two operands. If either operand is a /special value/--- then the general rules apply, except that if either operand is infinite and--- the other is finite an Invalid operation condition is raised and the result--- is NaN, or if both are infinite then the result is the first operand.------ Otherwise (both operands are finite), 'quantize' returns the number which--- is equal in value (except for any rounding) and sign to the first--- (left-hand) operand and which has an /exponent/ set to be equal to the--- exponent of the second (right-hand) operand.------ The /coefficient/ of the result is derived from that of the left-hand--- operand. It may be rounded using the current /rounding/ setting (if the--- /exponent/ is being increased), multiplied by a positive power of ten (if--- the /exponent/ is being decreased), or is unchanged (if the /exponent/ is--- already equal to that of the right-hand operand).------ Unlike other operations, if the length of the /coefficient/ after the--- quantize operation would be greater than /precision/ then an Invalid--- operation condition is raised. This guarantees that, unless there is an--- error condition, the /exponent/ of the result of a quantize is always equal--- to that of the right-hand operand.------ Also unlike other operations, quantize will never raise Underflow, even if--- the result is subnormal and inexact.-quantize :: (Precision p, Rounding r)-         => Decimal p r -> Decimal a b -> Arith p r (Decimal p r)-quantize x@Num { coefficient = xc, exponent = xe } Num { exponent = ye }-  | xe > ye   = result x { coefficient = xc * 10^(xe - ye), exponent = ye }-  | xe < ye   = rc >>= \c -> return x { coefficient = c, exponent = ye }-  | otherwise = return x--  where result :: Precision p => Decimal p r -> Arith p r (Decimal p r)-        result x = getPrecision >>= \p -> case numDigits (coefficient x) of-          n | maybe False (n >) p -> invalidOperation x-          _                       -> return x--        rc :: Rounding r => Arith p r Coefficient-        rc = let b      = 10^(ye - xe)-                 (q, r) = xc `quotRem` b-             in getRounder >>= \rounder -> return (rounder (sign x) r b q)--quantize Num{}      Inf{}    = invalidOperation qNaN-quantize Inf{}      Num{}    = invalidOperation qNaN-quantize n@Inf{}    Inf{}    = return n-quantize n@SNaN{}   _        = invalidOperation n-quantize _          n@SNaN{} = invalidOperation n-quantize n@QNaN{}   _        = return         n-quantize _          n@QNaN{} = return (coerce n)--{- $doctest-quantize->>> op2 Op.quantize "2.17" "0.001"-2.170-->>> op2 Op.quantize "2.17" "0.01"-2.17-->>> op2 Op.quantize "2.17" "0.1"-2.2-->>> op2 Op.quantize "2.17" "1e+0"-2-->>> op2 Op.quantize "2.17" "1e+1"-0E+1-->>> op2 Op.quantize "-Inf" "Infinity"--Infinity-->>> op2 Op.quantize "2" "Infinity"-NaN-->>> op2 Op.quantize "-0.1" "1"--0-->>> op2 Op.quantize "-0" "1e+5"--0E+5-->>> op2 Op.quantize "+35236450.6" "1e-2"-NaN-->>> op2 Op.quantize "-35236450.6" "1e-2"-NaN-->>> op2 Op.quantize "217" "1e-1"-217.0-->>> op2 Op.quantize "217" "1e+0"-217-->>> op2 Op.quantize "217" "1e+1"-2.2E+2-->>> op2 Op.quantize "217" "1e+2"-2E+2--}---- | 'reduce' takes one operand. It has the same semantics as the 'plus'--- operation, except that if the final result is finite it is reduced to its--- simplest form, with all trailing zeros removed and its sign preserved.-reduce :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r)-reduce n = reduce' <$> plus n-  where reduce' n@Num { coefficient = c, exponent = e }-          | c == 0 =         n {                  exponent = 0     }-          | r == 0 = reduce' n { coefficient = q, exponent = e + 1 }-          where (q, r) = c `quotRem` 10-        reduce' n = n--{- $doctest-reduce->>> op1 Op.reduce "2.1"-2.1-->>> op1 Op.reduce "-2.0"--2-->>> op1 Op.reduce "1.200"-1.2-->>> op1 Op.reduce "-120"--1.2E+2-->>> op1 Op.reduce "120.00"-1.2E+2-->>> op1 Op.reduce "0.00"-0--}---- | 'squareRoot' takes one operand. If the operand is a /special value/ then--- the general rules apply.------ Otherwise, the ideal exponent of the result is defined to be half the--- exponent of the operand (rounded to an integer, towards −Infinity, if--- necessary) and then:------ If the operand is less than zero an Invalid operation condition is raised.------ If the operand is greater than zero, the result is the square root of the--- operand. If no rounding is necessary (the exact result requires /precision/--- digits or fewer) then the coefficient and exponent giving the correct value--- and with the exponent closest to the ideal exponent is used. If the result--- must be inexact, it is rounded using the /round-half-even/ algorithm and--- the coefficient will have exactly /precision/ digits (unless the result is--- subnormal), and the exponent will be set to maintain the correct value.------ Otherwise (the operand is equal to zero), the result will be the zero with--- the same sign as the operand and with the ideal exponent.-squareRoot :: FinitePrecision p-           => Decimal a b -> Arith p r (Decimal p RoundHalfEven)-squareRoot n@Num { sign = s, coefficient = c, exponent = e }-  | c == 0   = return n { exponent = idealExp }-  | s == Pos = subResult >>= subRounded >>= result--  where idealExp = e `div` 2 :: Exponent--        reduced :: Decimal p r -> Decimal p r-        reduced n@Num { coefficient = c, exponent = e }-          | e < idealExp = case bd of-              Just (b, (q, _)) -> n { coefficient = q, exponent = e + b }-              Nothing          -> n-          | e > idealExp = n { coefficient = c * 10^d, exponent = idealExp }-          where d  = Prelude.abs (e - idealExp)-                bd = find (\(_, (_, r)) -> r == 0) ds-                ds = map (\d -> (d, c `quotRem` (10^d))) [d, d - 1 .. 1]-        reduced n = n--        subResult :: FinitePrecision p-                  => Arith p r (Decimal (PPlus1 (PPlus1 p)) RoundHalfEven)-        subResult = subArith (babylonian approx)--        subRounded :: Precision p-                   => Decimal a b -> Arith p r (Decimal p RoundHalfEven)-        subRounded = subArith . roundDecimal--        exactness :: Decimal a b -> Arith p r (Decimal PInfinite RoundHalfEven)-        exactness r = subArith (multiply r r >>= compare n)--        result :: Decimal p a -> Arith p r (Decimal p a)-        result r = do-          e <- exactness r-          if Number.isZero e-            then return (reduced r)-            else let r' = coerce r-                 in coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r')--        approx :: Decimal p r-        approx | even ae   = n { coefficient = 2, exponent =  ae      `quot` 2 }-               | otherwise = n { coefficient = 6, exponent = (ae - 1) `quot` 2 }-          where ae = adjustedExponent n--        babylonian :: FinitePrecision p => Decimal p RoundHalfEven-                   -> Arith p RoundHalfEven (Decimal p RoundHalfEven)-        babylonian x = do-          x' <- multiply oneHalf =<< add x =<< n `divide` x-          if x' == x then return x' else babylonian x'--squareRoot n@Inf { sign = Pos } = return (coerce n)-squareRoot n@QNaN{}             = return (coerce n)-squareRoot n                    = coerce <$> invalidOperation n--{- $doctest-squareRoot->>> op1 Op.squareRoot "0"-0-->>> op1 Op.squareRoot "-0"--0--This example appears to contradict the specification that the resulting-coefficient will have exactly /precision/ digits; awaiting clarification.-<<< op1 Op.squareRoot "0.39"-0.62449980-->>> op1 Op.squareRoot "100"-10-->>> op1 Op.squareRoot "1"-1-->>> op1 Op.squareRoot "1.0"-1.0-->>> op1 Op.squareRoot "1.00"-1.0-->>> op1 Op.squareRoot "7"-2.64575131-->>> op1 Op.squareRoot "10"-3.16227766--}---- $miscellaneous-operations------ This section describes miscellaneous operations on decimal numbers,--- including non-numeric comparisons, sign and other manipulations, and--- logical operations.------ The logical operations ('and', 'invert', 'or', and 'xor') take--- /logical operands/, which are finite numbers with a /sign/ of 0, an--- /exponent/ of 0, and a /coefficient/ whose digits must all be either 0 or--- 1. The length of the result will be at most /precision/ digits (all of--- which will be either 0 or 1); operands are truncated on the left or padded--- with zeros on the left as necessary. The result of a logical operation is--- never rounded and the only /flag/ that might be set is /invalid-operation/--- (set if an operand is not a valid logical operand).------ Some operations return a boolean value that is described as 0 or 1 in the--- documentation below. For reasons of efficiency, and as permitted by the--- /General Decimal Arithmetic Specification/, these operations return a--- 'Bool' in this implementation, but can be converted to 'Decimal' via--- 'fromBool'.--data Logical = Logical { bits :: Integer, bitLength :: Int }--toLogical :: Decimal a b -> Maybe Logical-toLogical Num { sign = Pos, coefficient = c, exponent = 0 } =-  getBits c Logical { bits = zeroBits, bitLength = 0 }--  where getBits :: Coefficient -> Logical -> Maybe Logical-        getBits 0 g = return g-        getBits c g@Logical { bits = b, bitLength = l } = case d of-          0 -> getBits c' g {                    bitLength = succ l }-          1 -> getBits c' g { bits = setBit b l, bitLength = succ l }-          _ -> Nothing-          where (c', d) = c `quotRem` 10--toLogical _ = Nothing--fromLogical :: Logical -> Decimal a b-fromLogical Logical { bits = b, bitLength = l } =-  Num { sign = Pos, coefficient = fromBits 0 1 0, exponent = 0 }--  where fromBits :: Int -> Coefficient -> Coefficient -> Coefficient-        fromBits i r c-          | i == l      = c-          | testBit b i = fromBits i' r' (c + r)-          | otherwise   = fromBits i' r'  c-          where i' = succ i-                r' = r * 10---- | 'and' is a logical operation which takes two logical operands. The result--- is the digit-wise /and/ of the two operands; each digit of the result is--- the logical and of the corresponding digits of the operands, aligned at the--- least-significant digit. A result digit is 1 if both of the corresponding--- operand digits are 1; otherwise it is 0.-and :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)-and x y = case (toLogical x, toLogical y) of-  (Just lx, Just ly) -> getPrecision >>= \p ->-    let m = Prelude.min (bitLength lx) (bitLength ly)-        z = Logical { bits = bits lx .&. bits ly-                    , bitLength = maybe m (Prelude.min m) p }-    in return (fromLogical z)-  _ -> invalidOperation qNaN--{- $doctest-and->>> op2 Op.and "0" "0"-0-->>> op2 Op.and "0" "1"-0-->>> op2 Op.and "1" "0"-0-->>> op2 Op.and "1" "1"-1-->>> op2 Op.and "1100" "1010"-1000-->>> op2 Op.and "1111" "10"-10--}---- | 'or' is a logical operation which takes two logical operands. The result--- is the digit-wise /inclusive or/ of the two operands; each digit of the--- result is the logical or of the corresponding digits of the operands,--- aligned at the least-significant digit. A result digit is 1 if either or--- both of the corresponding operand digits is 1; otherwise it is 0.-or :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)-or x y = case (toLogical x, toLogical y) of-  (Just lx, Just ly) -> getPrecision >>= \p ->-    let m = Prelude.max (bitLength lx) (bitLength ly)-        z = Logical { bits = bits lx .|. bits ly-                    , bitLength = maybe m (Prelude.min m) p }-    in return (fromLogical z)-  _ -> invalidOperation qNaN--{- $doctest-or->>> op2 Op.or "0" "0"-0-->>> op2 Op.or "0" "1"-1-->>> op2 Op.or "1" "0"-1-->>> op2 Op.or "1" "1"-1-->>> op2 Op.or "1100" "1010"-1110-->>> op2 Op.or "1110" "10"-1110--}---- | 'xor' is a logical operation which takes two logical operands. The result--- is the digit-wise /exclusive or/ of the two operands; each digit of the--- result is the logical exclusive-or of the corresponding digits of the--- operands, aligned at the least-significant digit. A result digit is 1 if--- one of the corresponding operand digits is 1 and the other is 0; otherwise--- it is 0.-xor :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)-xor x y = case (toLogical x, toLogical y) of-  (Just lx, Just ly) -> getPrecision >>= \p ->-    let m = Prelude.max (bitLength lx) (bitLength ly)-        z = Logical { bits = bits lx `Bits.xor` bits ly-                    , bitLength = maybe m (Prelude.min m) p }-    in return (fromLogical z)-  _ -> invalidOperation qNaN--{- $doctest-xor->>> op2 Op.xor "0" "0"-0-->>> op2 Op.xor "0" "1"-1-->>> op2 Op.xor "1" "0"-1-->>> op2 Op.xor "1" "1"-0-->>> op2 Op.xor "1100" "1010"-110-->>> op2 Op.xor "1111" "10"-1101--}---- | 'invert' is a logical operation which takes one logical operand. The--- result is the digit-wise /inversion/ of the operand; each digit of the--- result is the inverse of the corresponding digit of the operand. A result--- digit is 1 if the corresponding operand digit is 0; otherwise it is 0.-invert :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p r)-invert x = case toLogical x of-  Just lx -> getPrecision >>= \(Just p) ->-    let z = Logical { bits = complement (bits lx), bitLength = p }-    in return (fromLogical z)-  _ -> invalidOperation qNaN--{- $doctest-invert->>> op1 Op.invert "0"-111111111-->>> op1 Op.invert "1"-111111110-->>> op1 Op.invert "111111111"-0-->>> op1 Op.invert "101010101"-10101010--}---- | 'canonical' takes one operand. The result has the same value as the--- operand but always uses a /canonical/ encoding. The definition of--- /canonical/ is implementation-defined; if more than one internal encoding--- for a given NaN, Infinity, or finite number is possible then one--- “preferred” encoding is deemed canonical. This operation then returns the--- value using that preferred encoding.------ If all possible operands have just one internal encoding each, then--- 'canonical' always returns the operand unchanged (that is, it has the same--- effect as 'copy'). This operation is unaffected by context and is quiet —--- no /flags/ are changed in the context.-canonical :: Decimal a b -> Arith p r (Decimal a b)-canonical = return--{- $doctest-canonical->>> op1 Op.canonical "2.50"-2.50--}---- | 'class_' takes one operand. The result is an indication of the /class/ of--- the operand, where the class is one of ten possibilities, corresponding to--- one of the strings @"sNaN"@ (signaling NaN), @\"NaN"@ (quiet NaN),--- @"-Infinity"@ (negative infinity), @"-Normal"@ (negative normal finite--- number), @"-Subnormal"@ (negative subnormal finite number), @"-Zero"@--- (negative zero), @"+Zero"@ (non-negative zero), @"+Subnormal"@ (positive--- subnormal finite number), @"+Normal"@ (positive normal finite number), or--- @"+Infinity"@ (positive infinity). This operation is quiet; no /flags/ are--- changed in the context.------ Note that unlike the special values in the model, the sign of any NaN is--- ignored in the classification, as required by IEEE 754.-class_ :: Precision a => Decimal a b -> Arith p r Class-class_ n = return $ case n of-  Num {} | Number.isZero n      -> NumberClass (sign n) ZeroClass-         | Number.isSubnormal n -> NumberClass (sign n) SubnormalClass-         | otherwise            -> NumberClass (sign n) NormalClass-  Inf {}                        -> NumberClass (sign n) InfinityClass-  QNaN{}                        -> NaNClass QNaNClass-  SNaN{}                        -> NaNClass SNaNClass--data Class = NumberClass Sign NumberClass -- ^ Number (finite or infinite)-           | NaNClass NaNClass            -- ^ Not a number (quiet or signaling)-           deriving Eq--data NumberClass = ZeroClass       -- ^ Zero-                 | SubnormalClass  -- ^ Subnormal finite number-                 | NormalClass     -- ^ Normal finite number-                 | InfinityClass   -- ^ Infinity-                 deriving Eq--data NaNClass = QNaNClass  -- ^ Not a number (quiet)-              | SNaNClass  -- ^ Not a number (signaling)-              deriving Eq--instance Show Class where-  show c = case c of-    NumberClass s nc   -> signChar s : showNumberClass nc-    NaNClass QNaNClass ->       nan-    NaNClass SNaNClass -> 's' : nan--    where signChar :: Sign -> Char-          signChar Pos = '+'-          signChar Neg = '-'--          showNumberClass :: NumberClass -> String-          showNumberClass s = case s of-            ZeroClass      -> "Zero"-            SubnormalClass -> "Subnormal"-            NormalClass    -> "Normal"-            InfinityClass  -> "Infinity"--          nan :: String-          nan = "NaN"--{- $doctest-class_->>> op1 Op.class_ "Infinity"-+Infinity-->>> op1 Op.class_ "1E-10"-+Normal-->>> op1 Op.class_ "2.50"-+Normal-->>> op1 Op.class_ "0.1E-999"-+Subnormal-->>> op1 Op.class_ "0"-+Zero-->>> op1 Op.class_ "-0"--Zero-->>> op1 Op.class_ "-0.1E-999"--Subnormal-->>> op1 Op.class_ "-1E-10"--Normal-->>> op1 Op.class_ "-2.50"--Normal-->>> op1 Op.class_ "-Infinity"--Infinity-->>> op1 Op.class_ "NaN"-NaN-->>> op1 Op.class_ "-NaN"-NaN-->>> op1 Op.class_ "sNaN"-sNaN--}---- | 'copy' takes one operand. The result is a copy of the operand. This--- operation is unaffected by context and is quiet — no /flags/ are changed in--- the context.-copy :: Decimal a b -> Arith p r (Decimal a b)-copy = return--{- $doctest-copy->>> op1 Op.copy "2.1"-2.1-->>> op1 Op.copy "-1.00"--1.00--}---- | 'copyAbs' takes one operand. The result is a copy of the operand with the--- /sign/ set to 0. Unlike the 'abs' operation, this operation is unaffected--- by context and is quiet — no /flags/ are changed in the context.-copyAbs :: Decimal a b -> Arith p r (Decimal a b)-copyAbs n = return n { sign = Pos }--{- $doctest-copyAbs->>> op1 Op.copyAbs "2.1"-2.1-->>> op1 Op.copyAbs "-100"-100--}---- | 'copyNegate' takes one operand. The result is a copy of the operand with--- the /sign/ inverted (a /sign/ of 0 becomes 1 and vice versa). Unlike the--- 'minus' operation, this operation is unaffected by context and is quiet —--- no /flags/ are changed in the context.-copyNegate :: Decimal a b -> Arith p r (Decimal a b)-copyNegate n = return n { sign = negateSign (sign n) }--{- $doctest-copyNegate->>> op1 Op.copyNegate "101.5"--101.5-->>> op1 Op.copyNegate "-101.5"-101.5--}---- | 'copySign' takes two operands. The result is a copy of the first operand--- with the /sign/ set to be the same as the /sign/ of the second--- operand. This operation is unaffected by context and is quiet — no /flags/--- are changed in the context.-copySign :: Decimal a b -> Decimal c d -> Arith p r (Decimal a b)-copySign n m = return n { sign = sign m }--{- $doctest-copySign->>> op2 Op.copySign  "1.50"  "7.33"-1.50-->>> op2 Op.copySign "-1.50"  "7.33"-1.50-->>> op2 Op.copySign  "1.50" "-7.33"--1.50-->>> op2 Op.copySign "-1.50" "-7.33"--1.50--}---- | 'isCanonical' takes one operand. The result is 1 if the operand is--- /canonical/; otherwise it is 0. The definition of /canonical/ is--- implementation-defined; if more than one internal encoding for a given NaN,--- Infinity, or finite number is possible then one “preferred” encoding is--- deemed canonical. This operation then tests whether the internal encoding--- is that preferred encoding.------ If all possible operands have just one internal encoding each, then--- 'isCanonical' always returns 1. This operation is unaffected by context and--- is quiet — no /flags/ are changed in the context.-isCanonical :: Decimal a b -> Arith p r Bool-isCanonical _ = return True--{- $doctest-isCanonical->>> fromBool $ op1 Op.isCanonical "2.50"-1--}---- | 'isFinite' takes one operand. The result is 1 if the operand is neither--- infinite nor a NaN (that is, it is a normal number, a subnormal number, or--- a zero); otherwise it is 0. This operation is unaffected by context and is--- quiet — no /flags/ are changed in the context.-isFinite :: Decimal a b -> Arith p r Bool-isFinite = return . Number.isFinite--{- $doctest-isFinite->>> fromBool $ op1 Op.isFinite "2.50"-1-->>> fromBool $ op1 Op.isFinite "-0.3"-1-->>> fromBool $ op1 Op.isFinite "0"-1-->>> fromBool $ op1 Op.isFinite "Inf"-0-->>> fromBool $ op1 Op.isFinite "NaN"-0--}---- | 'isInfinite' takes one operand. The result is 1 if the operand is an--- Infinity; otherwise it is 0. This operation is unaffected by context and is--- quiet — no /flags/ are changed in the context.-isInfinite :: Decimal a b -> Arith p r Bool-isInfinite n = return $ case n of-  Inf{} -> True-  _     -> False--{- $doctest-isInfinite->>> fromBool $ op1 Op.isInfinite "2.50"-0-->>> fromBool $ op1 Op.isInfinite "-Inf"-1-->>> fromBool $ op1 Op.isInfinite "NaN"-0--}---- | 'isNaN' takes one operand. The result is 1 if the operand is a NaN (quiet--- or signaling); otherwise it is 0. This operation is unaffected by context--- and is quiet — no /flags/ are changed in the context.-isNaN :: Decimal a b -> Arith p r Bool-isNaN n = return $ case n of-  QNaN{} -> True-  SNaN{} -> True-  _      -> False--{- $doctest-isNaN->>> fromBool $ op1 Op.isNaN "2.50"-0-->>> fromBool $ op1 Op.isNaN "NaN"-1-->>> fromBool $ op1 Op.isNaN "-sNaN"-1--}---- | 'isNormal' takes one operand. The result is 1 if the operand is a--- positive or negative /normal number/; otherwise it is 0. This operation is--- quiet; no /flags/ are changed in the context.-isNormal :: Precision a => Decimal a b -> Arith p r Bool-isNormal = return . Number.isNormal--{- $doctest-isNormal->>> fromBool $ op1 Op.isNormal "2.50"-1-->>> fromBool $ op1 Op.isNormal "0.1E-999"-0-->>> fromBool $ op1 Op.isNormal "0.00"-0-->>> fromBool $ op1 Op.isNormal "-Inf"-0-->>> fromBool $ op1 Op.isNormal "NaN"-0--}---- | 'isQNaN' takes one operand. The result is 1 if the operand is a quiet--- NaN; otherwise it is 0. This operation is unaffected by context and is--- quiet — no /flags/ are changed in the context.-isQNaN :: Decimal a b -> Arith p r Bool-isQNaN n = return $ case n of-  QNaN{} -> True-  _      -> False--{- $doctest-isQNaN->>> fromBool $ op1 Op.isQNaN "2.50"-0-->>> fromBool $ op1 Op.isQNaN "NaN"-1-->>> fromBool $ op1 Op.isQNaN "sNaN"-0--}---- | 'isSigned' takes one operand. The result is 1 if the /sign/ of the--- operand is 1; otherwise it is 0. This operation is unaffected by context--- and is quiet — no /flags/ are changed in the context.-isSigned :: Decimal a b -> Arith p r Bool-isSigned = return . Number.isNegative--{- $doctest-isSigned->>> fromBool $ op1 Op.isSigned "2.50"-0-->>> fromBool $ op1 Op.isSigned "-12"-1-->>> fromBool $ op1 Op.isSigned "-0"-1--}---- | 'isSNaN' takes one operand. The result is 1 if the operand is a signaling--- NaN; otherwise it is 0. This operation is unaffected by context and is--- quiet — no /flags/ are changed in the context.-isSNaN :: Decimal a b -> Arith p r Bool-isSNaN n = return $ case n of-  SNaN{} -> True-  _      -> False--{- $doctest-isSNaN->>> fromBool $ op1 Op.isSNaN "2.50"-0-->>> fromBool $ op1 Op.isSNaN "NaN"-0-->>> fromBool $ op1 Op.isSNaN "sNaN"-1--}---- | 'isSubnormal' takes one operand. The result is 1 if the operand is a--- positive or negative /subnormal number/; otherwise it is 0. This operation--- is quiet; no /flags/ are changed in the context.-isSubnormal :: Precision a => Decimal a b -> Arith p r Bool-isSubnormal = return . Number.isSubnormal--{- $doctest-isSubnormal->>> fromBool $ op1 Op.isSubnormal "2.50"-0-->>> fromBool $ op1 Op.isSubnormal "0.1E-999"-1-->>> fromBool $ op1 Op.isSubnormal "0.00"-0-->>> fromBool $ op1 Op.isSubnormal "-Inf"-0-->>> fromBool $ op1 Op.isSubnormal "NaN"-0--}---- | 'isZero' takes one operand. The result is 1 if the operand is a zero;--- otherwise it is 0. This operation is unaffected by context and is quiet —--- no /flags/ are changed in the context.-isZero :: Decimal a b -> Arith p r Bool-isZero = return . Number.isZero--{- $doctest-isZero->>> fromBool $ op1 Op.isZero "0"-1-->>> fromBool $ op1 Op.isZero "2.50"-0-->>> fromBool $ op1 Op.isZero "-0E+2"-1--}---- | 'logb' takes one operand. If the operand is a NaN then the general--- arithmetic rules apply. If the operand is infinite then +Infinity is--- returned. If the operand is a zero, then −Infinity is returned and the--- Division by zero exceptional condition is raised.------ Otherwise, the result is the integer which is the exponent of the magnitude--- of the most significant digit of the operand (as though the operand were--- truncated to a single digit while maintaining the value of that digit and--- without limiting the resulting exponent). All results are exact unless an--- integer result does not fit in the available /precision/.-logb :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r)-logb Num { coefficient = c, exponent = e }-  | c == 0    = raiseSignal DivisionByZero Inf { sign = Neg }-  | otherwise = roundDecimal (fromInteger r :: Decimal PInfinite RoundHalfEven)-  where r = fromIntegral (numDigits c) - 1 + fromIntegral e :: Integer-logb Inf{} = return Inf { sign = Pos }-logb n@QNaN{} = return (coerce n)-logb n@SNaN{} = invalidOperation n--{- $doctest-logb->>> op1 Op.logb "250"-2-->>> op1 Op.logb "2.50"-0-->>> op1 Op.logb "0.03"--2-->>> op1 Op.logb "0"--Infinity--}---- | 'radix' takes no operands. The result is the radix (base) in which--- arithmetic is effected; for this specification the result will have the--- value 10.-radix :: Precision p => Arith p r (Decimal p r)-radix = return radix'-  where radix' = case precision radix' of-          Just 1 -> one { exponent    =  1 }-          _      -> one { coefficient = 10 }--{- $doctest-radix->>> op0 Op.radix-10--}---- | 'sameQuantum' takes two operands, and returns 1 if the two operands have--- the same /exponent/ or 0 otherwise. The result is never affected by either--- the sign or the coefficient of either operand.------ If either operand is a /special value/, 1 is returned only if both operands--- are NaNs or both are infinities.------ 'sameQuantum' does not change any /flags/ in the context.-sameQuantum :: Decimal a b -> Decimal c d -> Arith p r Bool-sameQuantum Num { exponent = e1 } Num { exponent = e2 }-  | e1 == e2  = return True-  | otherwise = return False-sameQuantum Inf {} Inf {} = return True-sameQuantum QNaN{} QNaN{} = return True-sameQuantum SNaN{} SNaN{} = return True-sameQuantum QNaN{} SNaN{} = return True-sameQuantum SNaN{} QNaN{} = return True-sameQuantum _      _      = return False--{- $doctest-sameQuantum->>> fromBool $ op2 Op.sameQuantum "2.17" "0.001"-0-->>> fromBool $ op2 Op.sameQuantum "2.17" "0.01"-1-->>> fromBool $ op2 Op.sameQuantum "2.17" "0.1"-0-->>> fromBool $ op2 Op.sameQuantum "2.17" "1"-0-->>> fromBool $ op2 Op.sameQuantum "Inf" "-Inf"-1-->>> fromBool $ op2 Op.sameQuantum "NaN" "NaN"-1--}---- | 'shift' takes two operands. The second operand must be an integer (with--- an /exponent/ of 0) in the range /−precision/ through /precision/. If the--- first operand is a NaN then the general arithmetic rules apply, and if it--- is infinite then the result is the Infinity unchanged.------ Otherwise (the first operand is finite) the result has the same /sign/ and--- /exponent/ as the first operand, and a /coefficient/ which is a shifted--- copy of the digits in the coefficient of the first operand. The number of--- places to shift is taken from the absolute value of the second operand,--- with the shift being to the left if the second operand is positive or to--- the right otherwise. Digits shifted into the coefficient are zeros.------ The only /flag/ that might be set is /invalid-operation/ (set if the first--- operand is an sNaN or the second is not valid).------ The 'rotate' operation can be used to rotate rather than shift a--- coefficient.-shift :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)-shift n@Num { coefficient = c } s@Num { sign = d, coefficient = sc }-  | validShift z s = return z-  where z = case precision z of-          Just p  -> y { coefficient = coefficient y `rem` 10 ^ p }-          Nothing -> y-        y = case d of-          Pos -> n { coefficient =  c  *     10 ^ sc }-          Neg -> n { coefficient =  c `quot` 10 ^ sc }-shift n@Inf{}  s | validShift z s = return z where z = coerce n-shift n@QNaN{} s | validShift z s = return z where z = coerce n-shift n        _                  = invalidOperation n--validShift :: Precision p => p -> Decimal a b -> Bool-validShift px Num { coefficient = c, exponent = 0 } =-  let p = fromIntegral <$> precision px in maybe True (c <=) p-validShift _ _ = False--{- $doctest-shift->>> op2 Op.shift "34" "8"-400000000-->>> op2 Op.shift "12" "9"-0-->>> op2 Op.shift "123456789" "-2"-1234567-->>> op2 Op.shift "123456789" "0"-123456789-->>> op2 Op.shift "123456789" "+2"-345678900--}---- | 'rotate' takes two operands. The second operand must be an integer (with--- an /exponent/ of 0) in the range /−precision/ through /precision/. If the--- first operand is a NaN then the general arithmetic rules apply, and if it--- is infinite then the result is the Infinity unchanged.------ Otherwise (the first operand is finite) the result has the same /sign/ and--- /exponent/ as the first operand, and a /coefficient/ which is a rotated--- copy of the digits in the coefficient of the first operand. The number of--- places of rotation is taken from the absolute value of the second operand,--- with the rotation being to the left if the second operand is positive or to--- the right otherwise.------ If the coefficient of the first operand has fewer than /precision/ digits,--- it is treated as though it were padded on the left with zeros to length--- /precision/ before the rotation. Similarly, if the coefficient of the first--- operand has more than /precision/ digits, it is truncated on the left--- before use.------ The only /flag/ that might be set is /invalid-operation/ (set if the first--- operand is an sNaN or the second is not valid).------ The 'shift' operation can be used to shift rather than rotate a--- coefficient.-rotate :: FinitePrecision p-       => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)-rotate n@Num { coefficient = c } s@Num { sign = d, coefficient = sc }-  | validShift z s = return z-  where z = n { coefficient = rc * b + (lc `rem` b) }-        (lc, rc) = c `quotRem` b'-        (b , b') = case d of-          Pos -> (10^sc , 10^sc')-          Neg -> (10^sc', 10^sc )-        Just p = precision z-        sc'    = p - fromIntegral sc-rotate n@Inf{}  s | validShift z s = return z where z = coerce n-rotate n@QNaN{} s | validShift z s = return z where z = coerce n-rotate n        _                  = invalidOperation n--{- $doctest-rotate->>> op2 Op.rotate "34" "8"-400000003-->>> op2 Op.rotate "12" "9"-12-->>> op2 Op.rotate "123456789" "-2"-891234567-->>> op2 Op.rotate "123456789" "0"-123456789-->>> op2 Op.rotate "123456789" "+2"-345678912--}+         -- ** General arithmetic+         add+       , subtract+       , multiply+       , divide+         -- divideInteger+         -- remainder+         -- remainderNear+       , power+       , squareRoot+       , fusedMultiplyAdd++         -- ** Exponential and logarithmic+       , exp+       , ln+       , log10++         -- ** Unary sign+       , plus+       , minus+       , abs++         -- ** Comparison+       , compare+       , compareSignal++       , min+       , max+       , minMagnitude+       , maxMagnitude++         -- ** Rounding and quantization++       , roundToIntegralValue+       , roundToIntegralExact+       , quantize+       , reduce++         -- nextMinus+         -- nextPlus+         -- nextToward++         -- * Miscellaneous operations+         -- $miscellaneous-operations++         -- ** Logic and shifting+         -- $logical-operations+       , and+       , or+       , xor+       , invert++       , shift+       , rotate++         -- ** Predicates+       , isZero+       , isSigned+       , isFinite+       , isInfinite+       , isNormal+       , isSubnormal+       , isNaN+       , isQNaN+       , isSNaN+       , isCanonical++         -- ** Total comparison and classification+       , compareTotal+       , compareTotalMagnitude++       , class', Class(..), Sign(..), NumberClass(..), NaNClass(..)++         -- ** Exponent manipulation+       , logb+       , scaleb+       , sameQuantum+       , radix++         -- ** Sign manipulation and conversion+       , copyAbs+       , copyNegate+       , copySign+       , copy++       , canonical+       ) where++import Prelude hiding (abs, and, compare, exp, exponent, isInfinite, isNaN,+                       max, min, or, subtract)+import qualified Prelude++import Control.Monad (join)+import Data.Bits (complement, setBit, testBit, zeroBits, (.&.), (.|.))+import Data.Coerce (coerce)+import Data.List (find)+import Data.Maybe (fromMaybe)++import qualified Data.Bits as Bits++import Numeric.Decimal.Arithmetic+import Numeric.Decimal.Exception+import Numeric.Decimal.Number hiding (isFinite, isNormal, isSubnormal, isZero)+import Numeric.Decimal.Precision+import Numeric.Decimal.Rounding++import qualified Numeric.Decimal.Number as Number++finitePrecision :: FinitePrecision p => Decimal p r -> Int+finitePrecision n = let Just p = precision n in p++roundingAlg :: Rounding r => Arith p r a -> RoundingAlgorithm+roundingAlg = rounding . arithRounding+  where arithRounding :: Arith p r a -> r+        arithRounding = undefined++result :: (Precision p, Rounding r) => Decimal p r -> Arith p r (Decimal p r)+result = roundDecimal  -- ...+--  | maybe False (numDigits c >) (precision r) = undefined++generalRules1 :: Decimal a b -> Arith p r (Decimal p r)+generalRules1 nan@NaN { signaling = False } = return (coerce nan)+generalRules1 x                             = invalidOperation x++generalRules2 :: Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+generalRules2 nan@NaN { signaling = True } _ = invalidOperation nan+generalRules2 _ nan@NaN { signaling = True } = invalidOperation nan+generalRules2 nan@NaN{} _                    = return (coerce nan)+generalRules2 _         nan@NaN{}            = return (coerce nan)+generalRules2 x         _                    = invalidOperation x++-- $arithmetic-operations+--+-- This section describes the arithmetic operations on, and some other+-- functions of, numbers, including subnormal numbers, negative zeros, and+-- special values (see also IEEE 754 §5 and §6).++-- | 'add' takes two operands. If either operand is a /special value/ then the+-- general rules apply.+--+-- Otherwise, the operands are added.+--+-- The result is then rounded to /precision/ digits if necessary, counting+-- from the most significant digit of the result.+add :: (Precision p, Rounding r)+    => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+add Num { sign = xs, coefficient = xc, exponent = xe }+    Num { sign = ys, coefficient = yc, exponent = ye } = sum++  where sum = result Num { sign = rs, coefficient = rc, exponent = re }+        rs | rc /= 0                       = if xac > yac then xs else ys+           | xs == Neg && ys == Neg        = Neg+           | xs /= ys &&+             roundingAlg sum == RoundFloor = Neg+           | otherwise                     = Pos+        rc | xs == ys  = xac + yac+           | xac > yac = xac - yac+           | otherwise = yac - xac+        re = Prelude.min xe ye+        (xac, yac) | xe == ye  = (xc, yc)+                   | xe >  ye  = (xc * 10^n, yc)+                   | otherwise = (xc, yc * 10^n)+          where n = Prelude.abs (xe - ye)++add inf@Inf{} Num{} = return (coerce inf)+add Num{} inf@Inf{} = return (coerce inf)+add inf@Inf { sign = xs } Inf { sign = ys }+  | xs == ys  = return (coerce inf)+  | otherwise = invalidOperation qNaN+add x y = generalRules2 x y++-- | 'subtract' takes two operands. If either operand is a /special value/+-- then the general rules apply.+--+-- Otherwise, the operands are added after inverting the /sign/ used for the+-- second operand.+--+-- The result is then rounded to /precision/ digits if necessary, counting+-- from the most significant digit of the result.+subtract :: (Precision p, Rounding r)+         => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+subtract x@Num{} y@Num{} = add x (flipSign y)+subtract x@Inf{} y@Num{} = add x (flipSign y)+subtract x@Num{} y@Inf{} = add x (flipSign y)+subtract x@Inf{} y@Inf{} = add x (flipSign y)+subtract x y = generalRules2 x y++-- | 'minus' takes one operand, and corresponds to the prefix minus operator+-- in programming languages.+--+-- Note that the result of this operation is affected by context and may set+-- /flags/. The 'copyNegate' operation may be used instead of 'minus' if this+-- is not desired.+minus :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r)+minus x = zero { exponent = exponent x } `subtract` x++-- | 'plus' takes one operand, and corresponds to the prefix plus operator in+-- programming languages.+--+-- Note that the result of this operation is affected by context and may set+-- /flags/.+plus :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r)+plus x = zero { exponent = exponent x } `add` x++-- | 'multiply' takes two operands. If either operand is a /special value/+-- then the general rules apply. Otherwise, the operands are multiplied+-- together (“long multiplication”), resulting in a number which may be as+-- long as the sum of the lengths of the two operands.+--+-- The result is then rounded to /precision/ digits if necessary, counting+-- from the most significant digit of the result.+multiply :: (Precision p, Rounding r)+         => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+multiply Num { sign = xs, coefficient = xc, exponent = xe }+         Num { sign = ys, coefficient = yc, exponent = ye } = result rn++  where rn = Num { sign = rs, coefficient = rc, exponent = re }+        rs = xorSigns xs ys+        rc = xc * yc+        re = xe + ye++multiply Inf { sign = xs } Inf { sign = ys } =+  return Inf { sign = xorSigns xs ys }+multiply Inf { sign = xs } Num { sign = ys, coefficient = yc }+  | yc == 0   = invalidOperation qNaN+  | otherwise = return Inf { sign = xorSigns xs ys }+multiply Num { sign = xs, coefficient = xc } Inf { sign = ys }+  | xc == 0   = invalidOperation qNaN+  | otherwise = return Inf { sign = xorSigns xs ys }+multiply x y = generalRules2 x y++-- | 'exp' takes one operand. If the operand is a NaN then the general rules+-- for special values apply.+--+-- Otherwise, the result is /e/ raised to the power of the operand, with the+-- following cases:+--+-- * If the operand is −Infinity, the result is 0 and exact.+--+-- * If the operand is a zero, the result is 1 and exact.+--+-- * If the operand is +Infinity, the result is +Infinity and exact.+--+-- * Otherwise the result is inexact and will be rounded using the+-- 'RoundHalfEven' algorithm. The coefficient will have exactly /precision/+-- digits (unless the result is subnormal). These inexact results should be+-- correctly rounded, but may be up to 1 ulp (unit in last place) in error.+exp :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p RoundHalfEven)+exp x@Num { sign = s, coefficient = c }+  | c == 0    = return one+  | s == Neg  = subArith (maclaurin x { sign = Pos } >>= reciprocal) >>=+                subRounded >>= result+  | otherwise = subArith (maclaurin x) >>= subRounded >>= result++  where multiplyExact :: Decimal a b -> Decimal c d+                      -> Arith PInfinite RoundHalfEven+                         (Decimal PInfinite RoundHalfEven)+        multiplyExact = multiply++        maclaurin :: FinitePrecision p => Decimal a b+                  -> Arith p RoundHalfEven (Decimal p RoundHalfEven)+        maclaurin x+          | adjustedExponent x >= 0 = subArith (subMaclaurin x) >>= subRounded+          | otherwise = sum one one one one+          where sum :: FinitePrecision p+                    => Decimal p RoundHalfEven+                    -> Decimal PInfinite RoundHalfEven+                    -> Decimal PInfinite RoundHalfEven+                    -> Decimal PInfinite RoundHalfEven+                    -> Arith p RoundHalfEven (Decimal p RoundHalfEven)+                sum s num den n = do+                  num' <- subArith (multiplyExact num x)+                  den' <- subArith (multiplyExact den n)+                  s' <- add s =<< divide num' den'+                  if s' == s then return s'+                    else sum s' num' den' =<< subArith (add n one)++        subMaclaurin :: FinitePrecision p => Decimal a b+                     -> Arith p RoundHalfEven (Decimal p RoundHalfEven)+        subMaclaurin x = subArith (multiplyExact x oneHalf) >>= maclaurin >>=+          \r -> multiply r r++        subRounded :: Precision p+                   => Decimal (PPlus1 (PPlus1 p)) a+                   -> Arith p r (Decimal p RoundHalfEven)+        subRounded = subArith . roundDecimal++        result :: Decimal p a -> Arith p r (Decimal p a)+        result r = coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r')+          where r' = coerce r++exp n@Inf { sign = s }+  | s == Pos  = return (coerce n)+  | otherwise = return zero+exp x = coerce <$> generalRules1 x++-- | 'fusedMultiplyAdd' takes three operands; the first two are multiplied+-- together, using 'multiply', with sufficient precision and exponent range+-- that the result is exact and unrounded. No /flags/ are set by the+-- multiplication unless one of the first two operands is a signaling NaN or+-- one is a zero and the other is an infinity.+--+-- Unless the multiplication failed, the third operand is then added to the+-- result of that multiplication, using 'add', under the current context.+--+-- In other words, @fusedMultiplyAdd x y z@ delivers a result which is @(x ×+-- y) + z@ with only the one, final, rounding.+fusedMultiplyAdd :: (Precision p, Rounding r)+                 => Decimal a b -> Decimal c d -> Decimal e f+                 -> Arith p r (Decimal p r)+fusedMultiplyAdd x y z =+  either raise (return . coerce) (exactMult x y) >>= add z++  where exactMult :: Rounding r => Decimal a b -> Decimal c d+                  -> Either (Exception PInfinite r) (Decimal PInfinite r)+        exactMult x y = evalArith (multiply x y) newContext++        raise :: Exception a r -> Arith p r (Decimal p r)+        raise e = raiseSignal (exceptionSignal e) (coerce $ exceptionResult e)++-- | 'ln' takes one operand. If the operand is a NaN then the general rules+-- for special values apply.+--+-- Otherwise, the operand must be a zero or positive, and the result is the+-- natural (base /e/) logarithm of the operand, with the following cases:+--+-- * If the operand is a zero, the result is −Infinity and exact.+--+-- * If the operand is +Infinity, the result is +Infinity and exact.+--+-- * If the operand equals one, the result is 0 and exact.+--+-- * Otherwise the result is inexact and will be rounded using the+-- 'RoundHalfEven' algorithm. The coefficient will have exactly /precision/+-- digits (unless the result is subnormal). These inexact results should be+-- correctly rounded, but may be up to 1 ulp (unit in last place) in error.+ln :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p RoundHalfEven)+ln x@Num { sign = s, coefficient = c, exponent = e }+  | c == 0   = return infinity { sign = Neg }+  | s == Pos = if e <= 0 && c == 10^(-e) then return zero+               else subArith (subLn x) >>= subRounded >>= result++  where subLn :: FinitePrecision p => Decimal a b+              -> Arith p RoundHalfEven (Decimal p RoundHalfEven)+        subLn x = do+          let fe = fromIntegral (-(numDigits c - 1)) :: Exponent+              r  = fromIntegral (e - fe) :: Decimal PInfinite RoundHalfEven+          lnf <- taylorLn x { exponent = fe }+          add lnf =<< multiply r =<< ln10++        subRounded :: Precision p => Decimal (PPlus1 (PPlus1 p)) a+                   -> Arith p r (Decimal p RoundHalfEven)+        subRounded = subArith . roundDecimal++        result :: Decimal p a -> Arith p r (Decimal p a)+        result r = coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r')+          where r' = coerce r++ln n@Inf { sign = Pos } = return (coerce n)+ln x = coerce <$> generalRules1 x++taylorLn :: FinitePrecision p => Decimal a b+         -> Arith p RoundHalfEven (Decimal p RoundHalfEven)+taylorLn x = do+  num <- x `subtract` one+  den <- x `add`      one+  multiply two =<< sum =<< num `divide` den++    where sum :: FinitePrecision p => Decimal p RoundHalfEven+              -> Arith p RoundHalfEven (Decimal p RoundHalfEven)+          sum b = multiply b b >>= \b2 -> sum' b b b2 one++            where sum' :: FinitePrecision p+                       => Decimal p RoundHalfEven+                       -> Decimal p RoundHalfEven+                       -> Decimal p RoundHalfEven+                       -> Decimal PInfinite RoundHalfEven+                       -> Arith p RoundHalfEven (Decimal p RoundHalfEven)+                  sum' s m b n = do+                    m' <- multiply m b+                    n' <- subArith (add n two)+                    s' <- add s =<< divide m' n'+                    if s' == s then return s' else sum' s' m' b n'++ln10 :: FinitePrecision p => Arith p r (Decimal p RoundHalfEven)+ln10 = getPrecision >>= \(Just p) ->+  if p <= 50 then return fastLn10 else slowLn10++  where fastLn10 :: FinitePrecision p => Decimal p RoundHalfEven+        fastLn10 = 2.3025850929940456840179914546843642076011014886288++        slowLn10 :: FinitePrecision p => Arith p r (Decimal p RoundHalfEven)+        slowLn10 = subArith (taylorLn ten) >>= subRound++          where subRound :: Precision p => Decimal (PPlus1 (PPlus1 p)) a+                         -> Arith p r (Decimal p RoundHalfEven)+                subRound = subArith . roundDecimal++-- | 'log10' takes one operand. If the operand is a NaN then the general rules+-- for special values apply.+--+-- Otherwise, the operand must be a zero or positive, and the result is the+-- base 10 logarithm of the operand, with the following cases:+--+-- * If the operand is a zero, the result is −Infinity and exact.+--+-- * If the operand is +Infinity, the result is +Infinity and exact.+--+-- * If the operand equals an integral power of ten (including 10^0 and+-- negative powers) and there is sufficient /precision/ to hold the integral+-- part of the result, the result is an integer (with an exponent of 0) and+-- exact.+--+-- * Otherwise the result is inexact and will be rounded using the+-- 'RoundHalfEven' algorithm. The coefficient will have exactly /precision/+-- digits (unless the result is subnormal). These inexact results should be+-- correctly rounded, but may be up to 1 ulp (unit in last place) in error.+log10 :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p RoundHalfEven)+log10 x@Num { sign = s, coefficient = c, exponent = e }+  | c == 0   = return infinity { sign = Neg }+  | s == Pos = getPrecision >>= \prec -> case powerOfTen c e of+      Just p | maybe True (numDigits pc <=) prec -> return (fromInteger p)+        where pc = fromInteger (Prelude.abs p) :: Coefficient+      _ -> subArith (join $ divide <$> ln x <*> ln10) >>= result++  where powerOfTen :: Coefficient -> Exponent -> Maybe Integer+        powerOfTen c e+          | c == 10^d = Just (fromIntegral e + fromIntegral d)+          | otherwise = Nothing+          where d = numDigits c - 1 :: Int++        result :: Decimal p a -> Arith p r (Decimal p a)+        result r = coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r')+          where r' = coerce r++log10 n@Inf { sign = Pos } = return (coerce n)+log10 x = coerce <$> generalRules1 x++-- | 'divide' takes two operands. If either operand is a /special value/ then+-- the general rules apply.+--+-- Otherwise, if the divisor is zero then either the Division undefined+-- condition is raised (if the dividend is zero) and the result is NaN, or the+-- Division by zero condition is raised and the result is an Infinity with a+-- sign which is the /exclusive or/ of the signs of the operands.+--+-- Otherwise, a “long division” is effected.+--+-- The result is then rounded to /precision/ digits, if necessary, according+-- to the /rounding/ algorithm and taking into account the remainder from the+-- division.+divide :: (FinitePrecision p, Rounding r)+       => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+divide dividend@Num{ sign = xs } Num { coefficient = 0, sign = ys }+  | Number.isZero dividend = divisionUndefined+  | otherwise              = divisionByZero infinity { sign = xorSigns xs ys }+divide Num { sign = xs, coefficient = xc, exponent = xe }+       Num { sign = ys, coefficient = yc, exponent = ye } = quotient++  where quotient = result =<< answer+        rn = Num { sign = rs, coefficient = rc, exponent = re }+        rs = xorSigns xs ys+        (rc, rem, dv, adjust) = longDivision xc yc (finitePrecision rn)+        re = xe - (ye + adjust)+        answer+          | rem == 0  = return rn+          | otherwise = roundDecimal $ case (rem * 2) `Prelude.compare` dv of+              LT -> rn { coefficient = rc * 10 + 1, exponent = re - 1 }+              EQ -> rn { coefficient = rc * 10 + 5, exponent = re - 1 }+              GT -> rn { coefficient = rc * 10 + 9, exponent = re - 1 }++divide Inf{} Inf{} = invalidOperation qNaN+divide Inf { sign = xs } Num { sign = ys } =+  return Inf { sign = xorSigns xs ys }+divide Num { sign = xs } Inf { sign = ys } =+  return zero { sign = xorSigns xs ys }+divide x y = generalRules2 x y++type Dividend  = Coefficient+type Divisor   = Coefficient+type Quotient  = Coefficient+type Remainder = Dividend++longDivision :: Dividend -> Divisor -> Int+             -> (Quotient, Remainder, Divisor, Exponent)+longDivision 0  dv _ = (0, 0, dv, 0)+longDivision dd dv p = step1 dd dv 0++  where step1 :: Dividend -> Divisor -> Exponent+              -> (Quotient, Remainder, Divisor, Exponent)+        step1 dd dv adjust+          | dd <       dv = step1 (dd * 10)  dv       (adjust + 1)+          | dd >= 10 * dv = step1  dd       (dv * 10) (adjust - 1)+          | otherwise     = step2  dd        dv        adjust++        step2 :: Dividend -> Divisor -> Exponent+              -> (Quotient, Remainder, Divisor, Exponent)+        step2 = step3 0++        step3 :: Quotient -> Dividend -> Divisor -> Exponent+              -> (Quotient, Remainder, Divisor, Exponent)+        step3 r dd dv adjust+          | dv <= dd                 = step3 (r +  1) (dd - dv) dv  adjust+          | (dd == 0 && adjust >= 0) ||+            numDigits r == p         = step4  r        dd       dv  adjust+          | otherwise                = step3 (r * 10) (dd * 10) dv (adjust + 1)++        step4 :: Quotient -> Remainder -> Divisor -> Exponent+              -> (Quotient, Remainder, Divisor, Exponent)+        step4 = (,,,)++reciprocal :: (FinitePrecision p, Rounding r)+           => Decimal a b -> Arith p r (Decimal p r)+reciprocal = divide one++-- | 'abs' takes one operand. If the operand is negative, the result is the+-- same as using the 'minus' operation on the operand. Otherwise, the result+-- is the same as using the 'plus' operation on the operand.+--+-- Note that the result of this operation is affected by context and may set+-- /flags/. The 'copyAbs' operation may be used if this is not desired.+abs :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r)+abs x+  | isNegative x = minus x+  | otherwise    = plus  x++-- | 'compare' takes two operands and compares their values numerically. If+-- either operand is a /special value/ then the general rules apply. No flags+-- are set unless an operand is a signaling NaN.+--+-- Otherwise, the operands are compared, returning @'Right' 'LT'@ if the first+-- is less than the second, @'Right' 'EQ'@ if they are equal, or @'Right'+-- 'GT'@ if the first is greater than the second.+--+-- A 'Left' value is returned if the result is NaN, indicating an “unordered”+-- comparison (see IEEE 754 §5.11).+compare :: Decimal a b -> Decimal c d+        -> Arith p r (Either (Decimal p r) Ordering)+compare x@Num{} y@Num{} = nzp <$> subArith (subtract' xn yn)++  where subtract' :: Decimal a b -> Decimal c d+                  -> Arith PInfinite RoundHalfEven+                     (Decimal PInfinite RoundHalfEven)+        subtract' = subtract++        (xn, yn) | sign x /= sign y = (either id fromOrdering $ nzp x,+                                       either id fromOrdering $ nzp y)+                 | otherwise        = (x, y)++        nzp :: Decimal a b -> Either (Decimal p r) Ordering+        nzp Num { sign = s, coefficient = c }+          | c == 0    = Right EQ+          | s == Pos  = Right GT+          | otherwise = Right LT+        nzp Inf { sign = s } = case s of+          Pos -> Right GT+          Neg -> Right LT+        nzp n = Left (coerce n)++compare Inf { sign = xs } Inf { sign = ys } = return $ case (xs, ys) of+  (Pos, Neg) -> Right GT+  (Neg, Pos) -> Right LT+  _          -> Right EQ+compare Inf { sign = xs } Num { } = return $ case xs of+  Pos -> Right GT+  Neg -> Right LT+compare Num { } Inf { sign = ys } = return $ case ys of+  Pos -> Right LT+  Neg -> Right GT+compare x y = Left <$> generalRules2 x y++-- | 'compareSignal' takes two operands and compares their values+-- numerically. This operation is identical to 'compare', except that if+-- neither operand is a signaling NaN then any quiet NaN operand is treated as+-- though it were a signaling NaN. (That is, all NaNs signal, with signaling+-- NaNs taking precedence over quiet NaNs.)+compareSignal :: Decimal a b -> Decimal c d+              -> Arith p r (Either (Decimal p r) Ordering)+compareSignal x@NaN { signaling = True } y =     x `compare`     y+compareSignal x y@NaN { signaling = True } =     x `compare`     y+compareSignal x y                          = q2s x `compare` q2s y++  where q2s :: Decimal p r -> Decimal p r+        q2s nan@NaN{} = nan { signaling = True }+        q2s x         = x++-- | 'max' takes two operands, compares their values numerically, and returns+-- the maximum. If either operand is a NaN then the general rules apply,+-- unless one is a quiet NaN and the other is numeric, in which case the+-- numeric operand is returned.+max :: Decimal a b -> Decimal a b -> Arith p r (Decimal a b)+max x y = snd <$> minMax id x y++-- | 'maxMagnitude' takes two operands and compares their values numerically+-- with their /sign/ ignored and assumed to be 0.+--+-- If, without signs, the first operand is the larger then the original first+-- operand is returned (that is, with the original sign). If, without signs,+-- the second operand is the larger then the original second operand is+-- returned. Otherwise the result is the same as from the 'max' operation.+maxMagnitude :: Decimal a b -> Decimal a b -> Arith p r (Decimal a b)+maxMagnitude x y = snd <$> minMax withoutSign x y++-- | 'min' takes two operands, compares their values numerically, and returns+-- the minimum. If either operand is a NaN then the general rules apply,+-- unless one is a quiet NaN and the other is numeric, in which case the+-- numeric operand is returned.+min :: Decimal a b -> Decimal a b -> Arith p r (Decimal a b)+min x y = fst <$> minMax id x y++-- | 'minMagnitude' takes two operands and compares their values numerically+-- with their /sign/ ignored and assumed to be 0.+--+-- If, without signs, the first operand is the smaller then the original first+-- operand is returned (that is, with the original sign). If, without signs,+-- the second operand is the smaller then the original second operand is+-- returned. Otherwise the result is the same as from the 'min' operation.+minMagnitude :: Decimal a b -> Decimal a b -> Arith p r (Decimal a b)+minMagnitude x y = fst <$> minMax withoutSign x y++-- | Ordering function for 'min', 'minMagnitude', 'max', and 'maxMagnitude':+-- returns the original arguments as (smaller, larger) when the given function+-- is applied to them.+minMax :: (Decimal a b -> Decimal a b) -> Decimal a b -> Decimal a b+       -> Arith p r (Decimal a b, Decimal a b)+minMax _ x@Num{}                       NaN { signaling = False } = return (x, x)+minMax _ x@Inf{}                       NaN { signaling = False } = return (x, x)+minMax _   NaN { signaling = False } y@Num{}                     = return (y, y)+minMax _   NaN { signaling = False } y@Inf{}                     = return (y, y)++minMax f x y = f x `compare` f y >>= \c -> return $ case c of+  Right LT -> (x, y)+  Right GT -> (y, x)+  Right EQ -> case (sign x, sign y) of+    (Neg, Pos) -> (x, y)+    (Pos, Neg) -> (y, x)+    (Pos, Pos) -> case (x, y) of+      (Num { exponent = xe }, Num { exponent = ye }) | xe > ye -> (y, x)+      _ -> (x, y)+    (Neg, Neg) -> case (x, y) of+      (Num { exponent = xe }, Num { exponent = ye }) | xe < ye -> (y, x)+      _ -> (x, y)+  Left nan -> let nan' = coerce nan in (nan', nan')++withoutSign :: Decimal p r -> Decimal p r+withoutSign n = n { sign = Pos }++-- | 'power' takes two operands, and raises a number (the left-hand operand)+-- to a power (the right-hand operand). If either operand is a /special value/+-- then the general rules apply, except as stated below.+--+-- The following rules apply:+--+-- * If both operands are zero, or if the left-hand operand is less than zero+-- and the right-hand operand does not have an integral value or is infinite,+-- an Invalid operation condition is raised, the result is NaN, and the+-- following rules do not apply.+--+-- * If the left-hand operand is infinite, the result will be exact and will+-- be infinite if the right-hand side is positive, 1 if the right-hand side is+-- a zero, and 0 if the right-hand side is negative.+--+-- * If the left-hand operand is a zero, the result will be exact and will be+-- infinite if the right-hand side is negative or 0 if the right-hand side is+-- positive.+--+-- * If the right-hand operand is a zero, the result will be 1 and exact.+--+-- * In cases not covered above, the result will be inexact unless the+-- right-hand side has an integral value and the result is finite and can be+-- expressed exactly within /precision/ digits. In this latter case, if the+-- result is unrounded then its exponent will be that which would result if+-- the operation were calculated by repeated multiplication (if the second+-- operand is negative then the reciprocal of the first operand is used, with+-- the absolute value of the second operand determining the multiplications).+--+-- * Inexact finite results should be correctly rounded, but may be up to 1+-- ulp (unit in last place) in error.+--+-- * The /sign/ of the result will be 1 only if the right-hand side has an+-- integral value and is odd (and is not infinite) and also the /sign/ of the+-- left-hand side is 1. In all other cases, the /sign/ of the result will be+-- 0.+power :: (FinitePrecision p, Rounding r)+      => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+power x@Num { coefficient = 0 } y@Num{}+  | Number.isZero y     = invalidOperation qNaN+  | Number.isNegative y = divisionByZero infinity { sign = powerSign x y }+  | otherwise           = return zero     { sign = powerSign x y }+power x@Num{} y@Num{} = case integralValue y of+  Just i  | i < 0               -> reciprocal x >>= \rx -> integralPower rx (-i)+          | otherwise           ->                         integralPower  x   i+  Nothing | Number.isPositive x -> ln x >>= multiply y >>= fmap coerce . exp+          | otherwise           -> invalidOperation qNaN+power x@Num{} y@Inf{}+  | Number.isPositive x = return $ case sign y of+      Pos -> infinity+      Neg -> zero+  | otherwise           = invalidOperation qNaN+power x@Inf{} y@Num{}+  | Number.isZero y     = return one+  | Number.isPositive y = return infinity { sign = powerSign x y }+  | otherwise           = return zero     { sign = powerSign x y }+power Inf{} Inf { sign = s }+  | s == Pos            = return infinity+  | otherwise           = return zero+power x y = generalRules2 x y++powerSign :: Decimal a b -> Decimal c d -> Sign+powerSign x y+  | Number.isNegative x && fromMaybe False (odd <$> integralValue y) = Neg+  | otherwise                                                        = Pos++integralPower :: (Precision p, Rounding r)+              => Decimal a b -> Integer -> Arith p r (Decimal p r)+integralPower b e = integralPower' (return b) e one+  where integralPower' :: (Precision p, Rounding r)+                       => Arith p r (Decimal a b) -> Integer -> Decimal p r+                       -> Arith p r (Decimal p r)+        integralPower' _  0 r = return r+        integralPower' mb e r+          | odd e     = mb >>= \b -> multiply r b >>=+                        integralPower'              (multiply b b) e'+          | otherwise = integralPower' (mb >>= \b -> multiply b b) e' r+          where e' = e `div` 2++-- | 'quantize' takes two operands. If either operand is a /special value/+-- then the general rules apply, except that if either operand is infinite and+-- the other is finite an Invalid operation condition is raised and the result+-- is NaN, or if both are infinite then the result is the first operand.+--+-- Otherwise (both operands are finite), 'quantize' returns the number which+-- is equal in value (except for any rounding) and sign to the first+-- (left-hand) operand and which has an /exponent/ set to be equal to the+-- exponent of the second (right-hand) operand.+--+-- The /coefficient/ of the result is derived from that of the left-hand+-- operand. It may be rounded using the current /rounding/ setting (if the+-- /exponent/ is being increased), multiplied by a positive power of ten (if+-- the /exponent/ is being decreased), or is unchanged (if the /exponent/ is+-- already equal to that of the right-hand operand).+--+-- Unlike other operations, if the length of the /coefficient/ after the+-- quantize operation would be greater than /precision/ then an Invalid+-- operation condition is raised. This guarantees that, unless there is an+-- error condition, the /exponent/ of the result of a quantize is always equal+-- to that of the right-hand operand.+--+-- Also unlike other operations, quantize will never raise Underflow, even if+-- the result is subnormal and inexact.+quantize :: (Precision p, Rounding r)+         => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+quantize x@Num { coefficient = xc, exponent = xe } Num { exponent = ye }+  | xe > ye   = result x { coefficient = xc * 10^(xe - ye), exponent = ye }+  | xe < ye   = rc >>= \c -> result x { coefficient = c, exponent = ye }+  | otherwise = result x++  where result :: Precision p => Decimal a b -> Arith p r (Decimal p r)+        result x = getPrecision >>= \p -> case numDigits (coefficient x) of+          n | maybe False (n >) p -> invalidOperation qNaN+          _                       -> return (coerce x)++        rc :: Rounding r => Arith p r Coefficient+        rc = let b      = 10^(ye - xe)+                 (q, r) = xc `quotRem` b+             in getRounder >>= \rounder -> return (rounder (sign x) r b q)++quantize Num{}   Inf{} = invalidOperation qNaN+quantize Inf{}   Num{} = invalidOperation qNaN+quantize x@Inf{} Inf{} = return (coerce x)+quantize x y = generalRules2 x y++-- | 'reduce' takes one operand. It has the same semantics as the 'plus'+-- operation, except that if the final result is finite it is reduced to its+-- simplest form, with all trailing zeros removed and its sign preserved.+reduce :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r)+reduce n = reduce' <$> plus n+  where reduce' n@Num { coefficient = c, exponent = e }+          | c == 0 =         n {                  exponent = 0     }+          | r == 0 = reduce' n { coefficient = q, exponent = e + 1 }+          where (q, r) = c `quotRem` 10+        reduce' n = n++-- | 'roundToIntegralExact' takes one operand. If the operand is a+-- /special value/, or the exponent of the operand is non-negative, then the+-- result is the same as the operand (unless the operand is a signaling NaN,+-- as usual).+--+-- Otherwise (the operand has a negative exponent) the result is the same as+-- using the 'quantize' operation using the given operand as the+-- left-hand-operand, 1E+0 as the right-hand-operand, and the precision of the+-- operand as the /precision/ setting. The rounding mode is taken from the+-- context, as usual.+roundToIntegralExact :: (Precision a, Rounding r)+                     => Decimal a b -> Arith p r (Decimal a r)+roundToIntegralExact x@Num { exponent = e }+  | e >= 0    = return (coerce x)+  | otherwise =+      let (Right r, context) = runArith (quantize x one) newContext+          quantizeFlags = flags context++          maybeRaise :: Signal -> Decimal a r -> Arith p r (Decimal a r)+          maybeRaise sig+            | sig `signalMember` quantizeFlags =+                fmap coerce . raiseSignal sig . coerce+            | otherwise = return++      in maybeRaise Inexact r >> maybeRaise Rounded r++roundToIntegralExact x@Inf{} = return (coerce x)+roundToIntegralExact x = coerce <$> generalRules1 x++-- | 'roundToIntegralValue' takes one operand. It is identical to the+-- 'roundToIntegralExact' operation except that the 'Inexact' and 'Rounded'+-- flags are never set even if the operand is rounded (that is, the operation+-- is quiet unless the operand is a signaling NaN).+roundToIntegralValue :: (Precision a, Rounding r)+                     => Decimal a b -> Arith p r (Decimal a r)+roundToIntegralValue x@Num { exponent = e }+  | e >= 0    = return (coerce x)+  | otherwise = subArith (quantize x one)+roundToIntegralValue x@Inf{} = return (coerce x)+roundToIntegralValue x = coerce <$> generalRules1 x++-- | 'squareRoot' takes one operand. If the operand is a /special value/ then+-- the general rules apply.+--+-- Otherwise, the ideal exponent of the result is defined to be half the+-- exponent of the operand (rounded to an integer, towards −Infinity, if+-- necessary) and then:+--+-- If the operand is less than zero an Invalid operation condition is raised.+--+-- If the operand is greater than zero, the result is the square root of the+-- operand. If no rounding is necessary (the exact result requires /precision/+-- digits or fewer) then the coefficient and exponent giving the correct value+-- and with the exponent closest to the ideal exponent is used. If the result+-- must be inexact, it is rounded using the 'RoundHalfEven' algorithm and the+-- coefficient will have exactly /precision/ digits (unless the result is+-- subnormal), and the exponent will be set to maintain the correct value.+--+-- Otherwise (the operand is equal to zero), the result will be the zero with+-- the same sign as the operand and with the ideal exponent.+squareRoot :: FinitePrecision p+           => Decimal a b -> Arith p r (Decimal p RoundHalfEven)+squareRoot n@Num { sign = s, coefficient = c, exponent = e }+  | c == 0   = return n { exponent = idealExp }+  | s == Pos = subResult >>= subRounded >>= result++  where idealExp = e `div` 2 :: Exponent++        reduced :: Decimal p r -> Decimal p r+        reduced n@Num { coefficient = c, exponent = e }+          | e < idealExp = case bd of+              Just (b, (q, _)) -> n { coefficient = q, exponent = e + b }+              Nothing          -> n+          | e > idealExp = n { coefficient = c * 10^d, exponent = idealExp }+          where d  = Prelude.abs (e - idealExp)+                bd = find (\(_, (_, r)) -> r == 0) ds+                ds = map (\d -> (d, c `quotRem` (10^d))) [d, d - 1 .. 1]+        reduced n = n++        subResult :: FinitePrecision p+                  => Arith p r (Decimal (PPlus1 (PPlus1 p)) RoundHalfEven)+        subResult = subArith (babylonian approx)++        subRounded :: Precision p+                   => Decimal a b -> Arith p r (Decimal p RoundHalfEven)+        subRounded = subArith . roundDecimal++        exactness :: Decimal a b -> Arith p r+                     (Either (Decimal p r) Ordering)+        exactness r = subArith (multiply' r r) >>= compare n+          where multiply' :: Decimal a b -> Decimal c d+                          -> Arith PInfinite RoundHalfEven+                             (Decimal PInfinite RoundHalfEven)+                multiply' = multiply++        result :: Decimal p a -> Arith p r (Decimal p a)+        result r = exactness r >>= \e -> case e of+          Right EQ -> return (reduced r)+          _ -> let r' = coerce r+               in coerce <$> (raiseSignal Rounded =<< raiseSignal Inexact r')++        approx :: Decimal p r+        approx | even ae   = n { coefficient = 2, exponent =  ae      `quot` 2 }+               | otherwise = n { coefficient = 6, exponent = (ae - 1) `quot` 2 }+          where ae = adjustedExponent n++        babylonian :: FinitePrecision p => Decimal p RoundHalfEven+                   -> Arith p RoundHalfEven (Decimal p RoundHalfEven)+        babylonian x = do+          x' <- multiply oneHalf =<< add x =<< n `divide` x+          if x' == x then return x' else babylonian x'++squareRoot n@Inf { sign = Pos } = return (coerce n)+squareRoot x = coerce <$> generalRules1 x++-- $miscellaneous-operations+--+-- This section describes miscellaneous operations on decimal numbers,+-- including non-numeric comparisons, sign and other manipulations, and+-- logical operations.+--+-- Some operations return a boolean value described as 0 or 1 in the+-- /General Decimal Arithmetic Specification/, but which is returned as a+-- 'Bool' in this implementation. These values can be converted to 'Decimal'+-- via 'fromBool'.+--+-- Similarly, the total ordering operations return an 'Ordering' value in this+-- implementation, but can be converted to 'Decimal' via 'fromOrdering'.++data Logical = Logical { bits :: Integer, bitLength :: Int }++toLogical :: Decimal a b -> Maybe Logical+toLogical Num { sign = Pos, coefficient = c, exponent = 0 } =+  getBits c Logical { bits = zeroBits, bitLength = 0 }++  where getBits :: Coefficient -> Logical -> Maybe Logical+        getBits 0 g = return g+        getBits c g@Logical { bits = b, bitLength = l } = case d of+          0 -> getBits c' g {                    bitLength = succ l }+          1 -> getBits c' g { bits = setBit b l, bitLength = succ l }+          _ -> Nothing+          where (c', d) = c `quotRem` 10++toLogical _ = Nothing++fromLogical :: Logical -> Decimal a b+fromLogical Logical { bits = b, bitLength = l } =+  Num { sign = Pos, coefficient = fromBits 0 1 0, exponent = 0 }++  where fromBits :: Int -> Coefficient -> Coefficient -> Coefficient+        fromBits i r c+          | i == l      = c+          | testBit b i = fromBits i' r' (c + r)+          | otherwise   = fromBits i' r'  c+          where i' = succ i+                r' = r * 10++-- $logical-operations+--+-- The logical operations ('and', 'or', 'xor', and 'invert') take+-- /logical operands/, which are finite numbers with a /sign/ of 0, an+-- /exponent/ of 0, and a /coefficient/ whose digits must all be either 0 or+-- 1. The length of the result will be at most /precision/ digits (all of+-- which will be either 0 or 1); operands are truncated on the left or padded+-- with zeros on the left as necessary. The result of a logical operation is+-- never rounded and the only /flag/ that might be set is 'InvalidOperation'+-- (set if an operand is not a valid logical operand).++-- | 'and' is a logical operation which takes two logical operands. The result+-- is the digit-wise /and/ of the two operands; each digit of the result is+-- the logical and of the corresponding digits of the operands, aligned at the+-- least-significant digit. A result digit is 1 if both of the corresponding+-- operand digits are 1; otherwise it is 0.+and :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+and x@Num{} y@Num{} = case (toLogical x, toLogical y) of+  (Just lx, Just ly) -> getPrecision >>= \p ->+    let m = Prelude.min (bitLength lx) (bitLength ly)+        z = Logical { bits = bits lx .&. bits ly+                    , bitLength = maybe m (Prelude.min m) p }+    in return (fromLogical z)+  _ -> invalidOperation qNaN+and x y = generalRules2 x y++-- | 'or' is a logical operation which takes two logical operands. The result+-- is the digit-wise /inclusive or/ of the two operands; each digit of the+-- result is the logical or of the corresponding digits of the operands,+-- aligned at the least-significant digit. A result digit is 1 if either or+-- both of the corresponding operand digits is 1; otherwise it is 0.+or :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+or x@Num{} y@Num{} = case (toLogical x, toLogical y) of+  (Just lx, Just ly) -> getPrecision >>= \p ->+    let m = Prelude.max (bitLength lx) (bitLength ly)+        z = Logical { bits = bits lx .|. bits ly+                    , bitLength = maybe m (Prelude.min m) p }+    in return (fromLogical z)+  _ -> invalidOperation qNaN+or x y = generalRules2 x y++-- | 'xor' is a logical operation which takes two logical operands. The result+-- is the digit-wise /exclusive or/ of the two operands; each digit of the+-- result is the logical exclusive-or of the corresponding digits of the+-- operands, aligned at the least-significant digit. A result digit is 1 if+-- one of the corresponding operand digits is 1 and the other is 0; otherwise+-- it is 0.+xor :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+xor x@Num{} y@Num{} = case (toLogical x, toLogical y) of+  (Just lx, Just ly) -> getPrecision >>= \p ->+    let m = Prelude.max (bitLength lx) (bitLength ly)+        z = Logical { bits = bits lx `Bits.xor` bits ly+                    , bitLength = maybe m (Prelude.min m) p }+    in return (fromLogical z)+  _ -> invalidOperation qNaN+xor x y = generalRules2 x y++-- | 'invert' is a logical operation which takes one logical operand. The+-- result is the digit-wise /inversion/ of the operand; each digit of the+-- result is the inverse of the corresponding digit of the operand. A result+-- digit is 1 if the corresponding operand digit is 0; otherwise it is 0.+invert :: FinitePrecision p => Decimal a b -> Arith p r (Decimal p r)+invert x@Num{} = case toLogical x of+  Just lx -> getPrecision >>= \(Just p) ->+    let z = Logical { bits = complement (bits lx), bitLength = p }+    in return (fromLogical z)+  _ -> invalidOperation qNaN+invert x = generalRules1 x++-- | 'canonical' takes one operand. The result has the same value as the+-- operand but always uses a /canonical/ encoding. The definition of+-- /canonical/ is implementation-defined; if more than one internal encoding+-- for a given NaN, Infinity, or finite number is possible then one+-- “preferred” encoding is deemed canonical. This operation then returns the+-- value using that preferred encoding.+--+-- If all possible operands have just one internal encoding each, then+-- 'canonical' always returns the operand unchanged (that is, it has the same+-- effect as 'copy'). This operation is unaffected by context and is quiet —+-- no /flags/ are changed in the context.+canonical :: Decimal a b -> Arith p r (Decimal a b)+canonical = copy++-- | 'class'' takes one operand. The result is an indication of the /class/ of+-- the operand, where the class is one of ten possibilities, corresponding to+-- one of the strings @"sNaN"@ (signaling NaN), @\"NaN"@ (quiet NaN),+-- @"-Infinity"@ (negative infinity), @"-Normal"@ (negative normal finite+-- number), @"-Subnormal"@ (negative subnormal finite number), @"-Zero"@+-- (negative zero), @"+Zero"@ (non-negative zero), @"+Subnormal"@ (positive+-- subnormal finite number), @"+Normal"@ (positive normal finite number), or+-- @"+Infinity"@ (positive infinity). This operation is quiet; no /flags/ are+-- changed in the context.+--+-- Note that unlike the special values in the model, the sign of any NaN is+-- ignored in the classification, as required by IEEE 754.+class' :: Precision a => Decimal a b -> Arith p r Class+class' n = return $ case n of+  Num {} | Number.isZero n      -> NumberClass (sign n) ZeroClass+         | Number.isSubnormal n -> NumberClass (sign n) SubnormalClass+         | otherwise            -> NumberClass (sign n) NormalClass+  Inf {}                        -> NumberClass (sign n) InfinityClass+  NaN { signaling = s }         -> NaNClass (toEnum . fromEnum $ s)++data Class = NumberClass Sign NumberClass -- ^ Number (finite or infinite)+           | NaNClass NaNClass            -- ^ Not a number (quiet or signaling)+           deriving Eq++data NumberClass = ZeroClass       -- ^ Zero+                 | SubnormalClass  -- ^ Subnormal finite number+                 | NormalClass     -- ^ Normal finite number+                 | InfinityClass   -- ^ Infinity+                 deriving Eq++data NaNClass = QuietClass      -- ^ Quiet NaN+              | SignalingClass  -- ^ Signaling NaN+              deriving (Eq, Enum)++instance Show Class where+  show c = case c of+    NumberClass s nc        -> signChar s : showNumberClass nc+    NaNClass QuietClass     ->       nan+    NaNClass SignalingClass -> 's' : nan++    where signChar :: Sign -> Char+          signChar Pos = '+'+          signChar Neg = '-'++          showNumberClass :: NumberClass -> String+          showNumberClass nc = case nc of+            ZeroClass      -> "Zero"+            SubnormalClass -> "Subnormal"+            NormalClass    -> "Normal"+            InfinityClass  -> "Infinity"++          nan :: String+          nan = "NaN"++-- | 'compareTotal' takes two operands and compares them using their abstract+-- representation rather than their numerical value. A /total ordering/ is+-- defined for all possible abstract representations, as described below. If+-- the first operand is lower in the total order than the second operand then+-- the result is 'LT', if the operands have the same abstract representation+-- then the result is 'EQ', and if the first operand is higher in the total+-- order than the second operand then the result is 'GT'. The total ordering+-- is defined as follows.+--+-- 1. The following items describe the ordering for representations whose+-- /sign/ is 0. If the /sign/ is 1, the order is reversed. A representation+-- with a /sign/ of 1 is always lower in the ordering than one with a /sign/+-- of 0.+--+-- 2. Numbers (representations which are not NaNs) are ordered such that a+-- larger numerical value is higher in the ordering. If two representations+-- have the same numerical value then the exponent is taken into account;+-- larger (more positive) exponents are higher in the ordering.+--+-- 3. All quiet NaNs are higher in the total ordering than all signaling NaNs.+--+-- 4. Quiet NaNs and signaling NaNs are ordered according to their /payload/;+-- a larger payload is higher in the ordering.+--+-- For example, the following values are ordered from lowest to highest: @-NaN+-- -sNaN -Infinity -127 -1 -1.00 -0 -0.000 0 1.2300 1.23 1E+9 Infinity sNaN+-- NaN NaN456@.+compareTotal :: Decimal a b -> Decimal c d -> Arith p r Ordering+compareTotal x y = return $ case (sign x, sign y) of+  (Pos, Pos) -> compareAbs x y+  (Neg, Neg) -> compareAbs y x+  (Neg, Pos) -> LT+  (Pos, Neg) -> GT++  where compareAbs :: Decimal a b -> Decimal c d -> Ordering+        compareAbs Num { coefficient = xc, exponent = xe }+                   Num { coefficient = yc, exponent = ye } =+          let (xac, yac) | xe == ye  = (xc, yc)+                         | xe >  ye  = (xc * 10^n, yc)+                         | otherwise = (xc, yc * 10^n)+              n = Prelude.abs (xe - ye)+          in Prelude.compare xac yac `mappend` Prelude.compare xe ye+        compareAbs Num{} Inf{} = LT+        compareAbs Inf{} Num{} = GT+        compareAbs Inf{} Inf{} = EQ+        compareAbs NaN { signaling = xs, payload = xp }+                   NaN { signaling = ys, payload = yp } =+          Prelude.compare ys xs `mappend` Prelude.compare xp yp+        compareAbs NaN{} _     = GT+        compareAbs _     NaN{} = LT++-- | 'compareTotalMagnitude' takes two operands and compares them using their+-- abstract representation rather than their numerical value and with their+-- /sign/ ignored and assumed to be 0. The result is identical to that+-- obtained by using 'compareTotal' on two operands which are the 'copyAbs'+-- copies of the operands to 'compareTotalMagnitude'.+compareTotalMagnitude :: Decimal a b -> Decimal c d -> Arith p r Ordering+compareTotalMagnitude x y = compareTotal x { sign = Pos } y { sign = Pos }++-- | 'copy' takes one operand. The result is a copy of the operand. This+-- operation is unaffected by context and is quiet — no /flags/ are changed in+-- the context.+copy :: Decimal a b -> Arith p r (Decimal a b)+copy = return++-- | 'copyAbs' takes one operand. The result is a copy of the operand with the+-- /sign/ set to 0. Unlike the 'abs' operation, this operation is unaffected+-- by context and is quiet — no /flags/ are changed in the context.+copyAbs :: Decimal a b -> Arith p r (Decimal a b)+copyAbs n = return n { sign = Pos }++-- | 'copyNegate' takes one operand. The result is a copy of the operand with+-- the /sign/ inverted (a /sign/ of 0 becomes 1 and vice versa). Unlike the+-- 'minus' operation, this operation is unaffected by context and is quiet —+-- no /flags/ are changed in the context.+copyNegate :: Decimal a b -> Arith p r (Decimal a b)+copyNegate n = return n { sign = negateSign (sign n) }++-- | 'copySign' takes two operands. The result is a copy of the first operand+-- with the /sign/ set to be the same as the /sign/ of the second+-- operand. This operation is unaffected by context and is quiet — no /flags/+-- are changed in the context.+copySign :: Decimal a b -> Decimal c d -> Arith p r (Decimal a b)+copySign n m = return n { sign = sign m }++-- | 'isCanonical' takes one operand. The result is 'True' if the operand is+-- /canonical/; otherwise it is 'False'. The definition of /canonical/ is+-- implementation-defined; if more than one internal encoding for a given NaN,+-- Infinity, or finite number is possible then one “preferred” encoding is+-- deemed canonical. This operation then tests whether the internal encoding+-- is that preferred encoding.+--+-- If all possible operands have just one internal encoding each, then+-- 'isCanonical' always returns 'True'. This operation is unaffected by+-- context and is quiet — no /flags/ are changed in the context.+isCanonical :: Decimal a b -> Arith p r Bool+isCanonical _ = return True++-- | 'isFinite' takes one operand. The result is 'True' if the operand is+-- neither infinite nor a NaN (that is, it is a normal number, a subnormal+-- number, or a zero); otherwise it is 'False'. This operation is unaffected+-- by context and is quiet — no /flags/ are changed in the context.+isFinite :: Decimal a b -> Arith p r Bool+isFinite = return . Number.isFinite++-- | 'isInfinite' takes one operand. The result is 'True' if the operand is an+-- Infinity; otherwise it is 'False'. This operation is unaffected by context+-- and is quiet — no /flags/ are changed in the context.+isInfinite :: Decimal a b -> Arith p r Bool+isInfinite n = return $ case n of+  Inf{} -> True+  _     -> False++-- | 'isNaN' takes one operand. The result is 'True' if the operand is a NaN+-- (quiet or signaling); otherwise it is 'False'. This operation is unaffected+-- by context and is quiet — no /flags/ are changed in the context.+isNaN :: Decimal a b -> Arith p r Bool+isNaN n = return $ case n of+  NaN{} -> True+  _     -> False++-- | 'isNormal' takes one operand. The result is 'True' if the operand is a+-- positive or negative /normal number/; otherwise it is 'False'. This+-- operation is quiet; no /flags/ are changed in the context.+isNormal :: Precision a => Decimal a b -> Arith p r Bool+isNormal = return . Number.isNormal++-- | 'isQNaN' takes one operand. The result is 'True' if the operand is a+-- quiet NaN; otherwise it is 'False'. This operation is unaffected by context+-- and is quiet — no /flags/ are changed in the context.+isQNaN :: Decimal a b -> Arith p r Bool+isQNaN n = return $ case n of+  NaN { signaling = False } -> True+  _                         -> False++-- | 'isSigned' takes one operand. The result is 'True' if the /sign/ of the+-- operand is 1; otherwise it is 'False'. This operation is unaffected by+-- context and is quiet — no /flags/ are changed in the context.+isSigned :: Decimal a b -> Arith p r Bool+isSigned = return . Number.isNegative++-- | 'isSNaN' takes one operand. The result is 'True' if the operand is a+-- signaling NaN; otherwise it is 'False'. This operation is unaffected by+-- context and is quiet — no /flags/ are changed in the context.+isSNaN :: Decimal a b -> Arith p r Bool+isSNaN n = return $ case n of+  NaN { signaling = True } -> True+  _                        -> False++-- | 'isSubnormal' takes one operand. The result is 'True' if the operand is a+-- positive or negative /subnormal number/; otherwise it is 'False'. This+-- operation is quiet; no /flags/ are changed in the context.+isSubnormal :: Precision a => Decimal a b -> Arith p r Bool+isSubnormal = return . Number.isSubnormal++-- | 'isZero' takes one operand. The result is 'True' if the operand is a+-- zero; otherwise it is 'False'. This operation is unaffected by context and+-- is quiet — no /flags/ are changed in the context.+isZero :: Decimal a b -> Arith p r Bool+isZero = return . Number.isZero++-- | 'logb' takes one operand. If the operand is a NaN then the general+-- arithmetic rules apply. If the operand is infinite then +Infinity is+-- returned. If the operand is a zero, then −Infinity is returned and the+-- Division by zero exceptional condition is raised.+--+-- Otherwise, the result is the integer which is the exponent of the magnitude+-- of the most significant digit of the operand (as though the operand were+-- truncated to a single digit while maintaining the value of that digit and+-- without limiting the resulting exponent). All results are exact unless an+-- integer result does not fit in the available /precision/.+logb :: (Precision p, Rounding r) => Decimal a b -> Arith p r (Decimal p r)+logb Num { coefficient = c, exponent = e }+  | c == 0    = raiseSignal DivisionByZero Inf { sign = Neg }+  | otherwise = roundDecimal (fromInteger r :: Decimal PInfinite RoundHalfEven)+  where r = fromIntegral (numDigits c) - 1 + fromIntegral e :: Integer+logb Inf{} = return Inf { sign = Pos }+logb x = generalRules1 x++-- | 'scaleb' takes two operands. If either operand is a NaN then the general+-- arithmetic rules apply. Otherwise, the second operand must be a finite+-- integer with an exponent of zero and in the range ±2 × (E/max/ ++-- /precision/) inclusive, where E/max/ is the largest value that can be+-- returned by the 'logb' operation at the same /precision/ setting. (If is is+-- not, the Invalid Operation condition is raised and the result is NaN.)+--+-- If the first operand is infinite then that Infinity is returned, otherwise+-- the result is the first operand modified by adding the value of the second+-- operand to its /exponent/. The result may Overflow or Underflow.+scaleb :: Decimal a b -> Decimal c d -> Arith p r (Decimal a b)+scaleb x@Num { exponent = e } s+  | validScale s = let Just i = integralValue s+                   in return x { exponent = e + fromInteger i }+                      -- XXX check for Overflow and Underflow+scaleb x@Inf{} s | validScale s = return x+scaleb x y = coerce <$> generalRules2 x y++validScale :: Decimal a b -> Bool+validScale Num { exponent = 0 } = True  -- XXX+validScale _                    = False++-- | 'radix' takes no operands. The result is the radix (base) in which+-- arithmetic is effected; for this specification the result will have the+-- value 10.+radix :: Precision p => Arith p r (Decimal p r)+radix = return radix'+  where radix' = case precision radix' of+          Just 1 -> one { exponent    =  1 }+          _      -> one { coefficient = 10 }++-- | 'sameQuantum' takes two operands, and returns 'True' if the two operands+-- have the same /exponent/ or 'False' otherwise. The result is never affected+-- by either the sign or the coefficient of either operand.+--+-- If either operand is a /special value/, 'True' is returned only if both+-- operands are NaNs or both are infinities.+--+-- 'sameQuantum' does not change any /flags/ in the context.+sameQuantum :: Decimal a b -> Decimal c d -> Arith p r Bool+sameQuantum Num { exponent = xe } Num { exponent = ye } = return (xe == ye)+sameQuantum Inf {               } Inf {               } = return True+sameQuantum NaN {               } NaN {               } = return True+sameQuantum _                     _                     = return False++-- | 'shift' takes two operands. The second operand must be an integer (with+-- an /exponent/ of 0) in the range /−precision/ through /precision/. If the+-- first operand is a NaN then the general arithmetic rules apply, and if it+-- is infinite then the result is the Infinity unchanged.+--+-- Otherwise (the first operand is finite) the result has the same /sign/ and+-- /exponent/ as the first operand, and a /coefficient/ which is a shifted+-- copy of the digits in the coefficient of the first operand. The number of+-- places to shift is taken from the absolute value of the second operand,+-- with the shift being to the left if the second operand is positive or to+-- the right otherwise. Digits shifted into the coefficient are zeros.+--+-- The only /flag/ that might be set is 'InvalidOperation' (set if the first+-- operand is an sNaN or the second is not valid).+--+-- The 'rotate' operation can be used to rotate rather than shift a+-- coefficient.+shift :: Precision p => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+shift n@Num { coefficient = c } s@Num { sign = d, coefficient = sc }+  | validShift z s = return z+  where z = case precision z of+          Just p  -> y { coefficient = coefficient y `rem` 10 ^ p }+          Nothing -> y+        y = case d of+          Pos -> n { coefficient =  c  *     10 ^ sc }+          Neg -> n { coefficient =  c `quot` 10 ^ sc }++shift n@Inf {                   } s | validShift z s = return z+  where z = coerce n+shift n@NaN { signaling = False } s | validShift z s = return z+  where z = coerce n+shift n@NaN { signaling = True  } _                  = invalidOperation n+shift _                           s                  = invalidOperation s++validShift :: Precision p => p -> Decimal a b -> Bool+validShift px Num { coefficient = c, exponent = 0 } =+  let p = fromIntegral <$> precision px in maybe True (c <=) p+validShift _ _ = False++-- | 'rotate' takes two operands. The second operand must be an integer (with+-- an /exponent/ of 0) in the range /−precision/ through /precision/. If the+-- first operand is a NaN then the general arithmetic rules apply, and if it+-- is infinite then the result is the Infinity unchanged.+--+-- Otherwise (the first operand is finite) the result has the same /sign/ and+-- /exponent/ as the first operand, and a /coefficient/ which is a rotated+-- copy of the digits in the coefficient of the first operand. The number of+-- places of rotation is taken from the absolute value of the second operand,+-- with the rotation being to the left if the second operand is positive or to+-- the right otherwise.+--+-- If the coefficient of the first operand has fewer than /precision/ digits,+-- it is treated as though it were padded on the left with zeros to length+-- /precision/ before the rotation. Similarly, if the coefficient of the first+-- operand has more than /precision/ digits, it is truncated on the left+-- before use.+--+-- The only /flag/ that might be set is 'InvalidOperation' (set if the first+-- operand is an sNaN or the second is not valid).+--+-- The 'shift' operation can be used to shift rather than rotate a+-- coefficient.+rotate :: FinitePrecision p+       => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+rotate n@Num { coefficient = c } s@Num { sign = d, coefficient = sc }+  | validShift z s = return z+  where z = n { coefficient = rc * b + (lc `rem` b) }+        (lc, rc) = c `quotRem` b'+        (b , b') = case d of+          Pos -> (10^sc , 10^sc')+          Neg -> (10^sc', 10^sc )+        sc' = finitePrecision z - fromIntegral sc++rotate n@Inf {                   } s | validShift z s = return z+  where z = coerce n+rotate n@NaN { signaling = False } s | validShift z s = return z+  where z = coerce n+rotate n@NaN { signaling = True  } _                  = invalidOperation n+rotate _                           s                  = invalidOperation s

@@ -11,6 +11,7 @@        , minus        , abs        , compare+       , compareTotal        , min        , max        , power@@ -48,12 +49,11 @@          => Decimal a b -> Arith p r (Decimal p r) abs      :: (Precision p, Rounding r)          => Decimal a b -> Arith p r (Decimal p r)-compare  :: (Precision p, Rounding r)-         => Decimal a b -> Decimal c d -> Arith p r (Decimal p r)-min      :: (Precision p, Rounding r)-         => Decimal a b -> Decimal a b -> Arith p r (Decimal a b)-max      :: (Precision p, Rounding r)-         => Decimal a b -> Decimal a b -> Arith p r (Decimal a b)+compare  :: Decimal a b -> Decimal c d+         -> Arith p r (Either (Decimal p r) Ordering)+compareTotal :: Decimal a b -> Decimal c d -> Arith p r Ordering+min      :: Decimal a b -> Decimal a b -> Arith p r (Decimal a b)+max      :: Decimal a b -> Decimal a b -> Arith p r (Decimal a b) power    :: (FinitePrecision p, Rounding r)          => Decimal a b -> Decimal c d -> Arith p r (Decimal p r) squareRoot :: FinitePrecision p

@@ -16,12 +16,27 @@        , PInfinite        ) where +import {-# SOURCE #-} Numeric.Decimal.Number+ -- | Precision indicates the maximum number of significant decimal digits a -- number may have. class Precision p where-  -- | Return the precision of the argument, or 'Nothing' if the precision is infinite.+  -- | Return the precision of the argument, or 'Nothing' if the precision is+  -- infinite.   precision :: p -> Maybe Int +  -- | Return the maximum exponent for a number in scientific notation with+  -- the given precision, or 'Nothing' if the exponent has no limit.+  eMax :: p -> Maybe Exponent+  eMax n = subtract 1 . (10 ^) . numDigits <$> base+    where mlength = precision n                       :: Maybe Int+          base    = (10 *) . fromIntegral <$> mlength :: Maybe Coefficient++  -- | Return the minimum exponent for a number in scientific notation with+  -- the given precision, or 'Nothing' if the exponent has no limit.+  eMin :: p -> Maybe Exponent+  eMin = fmap (1 -) . eMax+ -- | A subclass of precisions that are finite class Precision p => FinitePrecision p @@ -29,6 +44,7 @@ data PInfinite instance Precision PInfinite where   precision _ = Nothing+  eMax      _ = Nothing  -- | A precision of 1 significant digit data P1

@@ -1,7 +1,7 @@  module Numeric.Decimal.Rounding        ( RoundingAlgorithm(..)-       , Rounding(..)+       , Rounding(rounding)         , RoundDown        , RoundHalfUp@@ -14,6 +14,8 @@        , Round05Up         , getRounder+       , Rounder+        , roundDecimal        ) where @@ -21,9 +23,10 @@  import Data.Coerce (coerce) +import {-# SOURCE #-} Numeric.Decimal.Arithmetic+import {-# SOURCE #-} Numeric.Decimal.Exception import {-# SOURCE #-} Numeric.Decimal.Number import                Numeric.Decimal.Precision-import {-# SOURCE #-} Numeric.Decimal.Arithmetic  -- | A value representation of a rounding algorithm (cf. 'Rounding'). data RoundingAlgorithm = RoundDown@@ -53,7 +56,7 @@         getRounder' = return roundCoefficient  -- | Round a 'Decimal' to the precision of the arithmetic context using the--- rounding algorithm of the arithmetic context.+-- rounding mode of the arithmetic context. roundDecimal :: (Precision p, Rounding r)              => Decimal a b -> Arith p r (Decimal p r) roundDecimal n@Num { sign = s, coefficient = c, exponent = e } = do@@ -68,20 +71,23 @@           n'     = case excessDigits c' =<< p of             Nothing -> n { coefficient = c'          , exponent =      e' }             _       -> n { coefficient = c' `quot` 10, exponent = succ e' }-          rounded :: Decimal p r -> Arith p r (Decimal p r)-          rounded-            | r /= 0    = raiseSignal Inexact-            | otherwise = return-      raiseSignal Rounded =<< rounded n'  -- XXX check for overflow+      rounded =<< (if r /= 0 then inexact else return) n'+      -- XXX check for overflow      Nothing -> return (coerce n) -  where excessDigits :: Coefficient -> Int -> Maybe Int-        excessDigits c p | d > p     = Just (d - p)-                         | otherwise = Nothing-          where d = numDigits c :: Int+roundDecimal n@NaN { payload = p } = do+  prec <- getPrecision+  case excessDigits p =<< (pred <$> prec) of+    Just _  -> return n { payload = 0 }+    Nothing -> return (coerce n)  roundDecimal n = return (coerce n)++excessDigits :: Coefficient -> Int -> Maybe Int+excessDigits c p | d > p     = Just (d - p)+                 | otherwise = Nothing+  where d = numDigits c :: Int  -- Required algorithms... 

@@ -1,13 +1,13 @@ # This file was automatically generated by 'stack init'-# +# # Some commonly used options have been documented as comments in this file. # For advanced use and comprehensive documentation of the format, please see:-# http://docs.haskellstack.org/en/stable/yaml_configuration/+# https://docs.haskellstack.org/en/stable/yaml_configuration/  # Resolver to choose a 'specific' stackage snapshot or a compiler version. # A snapshot resolver dictates the compiler version and the set of packages # to be used for project dependencies. For example:-# +# # resolver: lts-3.5 # resolver: nightly-2015-09-21 # resolver: ghc-7.10.2@@ -15,11 +15,11 @@ # resolver: #  name: custom-snapshot #  location: "./custom-snapshot.yaml"-resolver: lts-6.9+resolver: lts-9.1  # User packages to be built. # Various formats can be used as shown in the example below.-# +# # packages: # - some-directory # - https://example.com/foo/bar/baz-0.0.2.tar.gz@@ -31,12 +31,12 @@ #  subdirs: #  - auto-update #  - wai-# +# # A package marked 'extra-dep: true' will only be built if demanded by a # non-dependency (i.e. a user package), and its test suites and benchmarks # will not be run. This is useful for tweaking upstream packages. packages:-- '.'+- . # Dependency packages to be pulled from upstream that are not in the resolver # (e.g., acme-missiles-0.3) extra-deps: []@@ -49,18 +49,18 @@  # Control whether we use the GHC we find on the path # system-ghc: true-# +# # Require a specific version of stack, using version ranges # require-stack-version: -any # Default-# require-stack-version: ">=1.1"-# +# require-stack-version: ">=1.5"+# # Override the architecture used by stack, especially useful on Windows # arch: i386 # arch: x86_64-# +# # Extra directories used by stack for building # extra-include-dirs: [/path/to/dir] # extra-lib-dirs: [/path/to/dir]-# +# # Allow a newer minor version of GHC than the snapshot specifies # compiler-check: newer-minor

@@ -0,0 +1,31 @@++module Arbitrary+       ( Arbitrary+       ) where++import Numeric.Decimal+import Test.QuickCheck++infinity :: (Precision p, Rounding r) => Decimal p r+infinity = read "Infinity"++instance (Precision p, Rounding r) => Arbitrary (Decimal p r) where+  arbitrary = frequency [(85, genNum), (10, genInf)]++genNum :: (Precision p, Rounding r) => Gen (Decimal p r)+genNum = do+  c <- choose (-(10^10), 10^10) :: Gen Integer+  e <- choose (-99, 99)         :: Gen Integer+  return $ read (show c ++ 'E' : show e)++genInf :: (Precision p, Rounding r) => Gen (Decimal p r)+genInf = do+  s <- elements [-1, 1]+  return (s * infinity)++genNaN :: (Precision p, Rounding r) => Gen (Decimal p r)+genNaN = oneof [nan "", nan "s"]+  where nan kind = do+          s <- elements ["", "-"]+          p <- choose (0, 10000) :: Gen Integer+          return $ read (s ++ kind ++ "NaN" ++ show p)

@@ -0,0 +1,25 @@++{-# LANGUAGE OverloadedStrings #-}++module Numeric.Decimal.EncodingSpec (spec) where++import Test.Hspec+import Data.Binary++import Numeric.Decimal.Encoding++spec :: Spec+spec = do+  it "encodes (-7.50) correctly" $+    encode (read "-7.50" :: Decimal64) `shouldBe`+    "\xA2\x30\x00\x00\x00\x00\x03\xD0"+  it "decodes (-7.50) correctly" $+    (decode "\xA2\x30\x00\x00\x00\x00\x03\xD0" :: Decimal64) `shouldBe` (-7.50)++  it "decodes Infinity correctly" $+    (decode "\x78\xFF\xFF\xFF\xFF\xFF\xFF\xFF" :: Decimal64) `shouldSatisfy`+    \x -> isInfinite x && signum x == 1++-- prop> decode (encode x) == (x :: Decimal32)+-- prop> decode (encode x) == (x :: Decimal64)+-- prop> decode (encode x) == (x :: Decimal128)

@@ -0,0 +1,165 @@++module Numeric.Decimal.NumberSpec (spec) where++import Arbitrary ()++import Test.Hspec+import Test.QuickCheck+import Numeric.Decimal++spec :: Spec+spec = do+  describe "Read" $ do+    it "reads with Just correctly (positive)" $+      (read "Just 123" :: Maybe GeneralDecimal) `shouldBe` Just 123+    it "reads with Just correctly (negative)" $+      (read "Just (-12.0)" :: Maybe GeneralDecimal) `shouldBe` Just (-12)++  describe "Ord" $ do+    it "satisfies (>) invariant" $+      property $ \x y ->+      x > y ==> max x y == x && max y x == (x :: BasicDecimal)+    it "satisfies (<) invariant" $+      property $ \x y ->+      x < y ==> min x y == x && min y x == (x :: BasicDecimal)++    it "satisfies `max` invariant with respect to 1st arg" $+      property $ \x y ->+      max x y == x ==> x >= (y :: BasicDecimal)+    it "satisfies `max` invariant with respect to 2nd arg" $+      property $ \x y ->+      max x y == y ==> y >= (x :: BasicDecimal)+    it "satisfies `min` invariant with respect to 1st arg" $+      property $ \x y ->+      min x y == x ==> x <= (y :: BasicDecimal)+    it "satisfies `min` invariant with respect to 2nd arg" $+      property $ \x y ->+      min x y == y ==> y <= (x :: BasicDecimal)++  describe "Enum" $ do+    it "enumerates `enumFromTo` precisely" $+      ([1.7 .. 5.7] :: [BasicDecimal]) `shouldBe` [1.7, 2.7, 3.7, 4.7, 5.7]++    it "enumerates `enumFromThenTo` precisely (ascending)" $+      ([0, 0.1 .. 2] :: [BasicDecimal]) `shouldBe`+      [  0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,+       1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0]+    it "enumerates `enumFromThenTo` precisely (descending)" $+      ([2, 1.9 .. 0] :: [BasicDecimal]) `shouldBe`+      [  2, 1.9, 1.8, 1.7, 1.6, 1.5, 1.4, 1.3, 1.2, 1.1,+       1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.0]++  describe "Num" $ do+    it "satisfies relation between (+) and (*)" $+      property $ \x ->+      x + x == x * (2 :: GeneralDecimal)+    it "satisfies relation between (-) and 0" $+      property $ \x ->+      isFinite x ==> x - x == (0 :: BasicDecimal)+    it "satisfies relation between (+) and `negate`" $+      property $ \x ->+      isFinite x ==> x + negate x == (0 :: BasicDecimal)+    it "satisfies relation between `abs` and 0" $+      property $ \x ->+      abs x >= (0 :: GeneralDecimal)+    it "satisfies relation between `abs` and `signum`" $+      property $ \x ->+      abs x * signum x == (x :: GeneralDecimal)++  describe "Fractional" $ do+    context "RoundHalfUp" $ do+      it "properly rounds (down)" $+        (4.14 :: Decimal P2 RoundHalfUp)   `shouldBe` 4.1+      it "properly rounds (up)" $+        (4.15 :: Decimal P2 RoundHalfUp)   `shouldBe` 4.2+    context "RoundHalfDown" $+      it "properly rounds" $+        (4.15 :: Decimal P2 RoundHalfDown) `shouldBe` 4.1+    context "RoundHalfEven" $ do+      it "properly rounds 1 (up)" $+        (4.15 :: Decimal P2 RoundHalfEven) `shouldBe` 4.2+      it "properly rounds 1 (down)" $+        (4.25 :: Decimal P2 RoundHalfEven) `shouldBe` 4.2+      it "properly rounds 2 (up)" $+        (4.35 :: Decimal P2 RoundHalfEven) `shouldBe` 4.4+      it "properly rounds 2 (down)" $+        (4.45 :: Decimal P2 RoundHalfEven) `shouldBe` 4.4++  describe "RealFrac" $ do+    it "satisfies `properFraction` invariant 1" $+      property $ \x ->+      let (n,f) = properFraction (x :: BasicDecimal) :: (Integer, BasicDecimal)+      in x == fromIntegral n + f+    it "satisfies `properFraction` invariant 2" $+      property $ \x ->+      let (n,_) = properFraction (x :: BasicDecimal) :: (Integer, BasicDecimal)+      in (x < 0 && n <= 0) || (x >= 0 && n >= 0)+    it "satisfies `properFraction` invariant 3" $+      property $ \x ->+      let (_,f) = properFraction (x :: BasicDecimal) :: (Integer, BasicDecimal)+      in (x < 0 && f <= 0) || (x >= 0 && f >= 0)+    it "satisfies `properFraction` invariant 4" $+      property $ \x ->+      let (_,f) = properFraction (x :: BasicDecimal) :: (Integer, BasicDecimal)+      in isFinite f ==> abs f < 1++  describe "Floating" $ do+    it "produces same `pi` as Double" $+      realToFrac (pi :: ExtendedDecimal P16) `shouldBe` (pi :: Double)+    it "satisfies relation between (**) and (^)" $+      property $ \x y ->+      y >= 0 ==> (x :: BasicDecimal) ** fromInteger y == x ^ y+    it "computes `sqrt` correctly" $+      property $ \x ->+      isFinite x && x >= 0 ==>+      sqrt (x * x) `shouldBe` (x :: ExtendedDecimal P32)+      -- coefficient (sqrt (x * x) - (x :: ExtendedDecimal P16)) <= 1++  describe "RealFloat" $ do+    it "satisfies `decodeFloat` invariant 1" $+      property $ \x ->+      isFinite x ==> let b      = floatRadix (x :: BasicDecimal)+                         (m, n) = decodeFloat x+                     in x == fromInteger m * fromInteger b ^^ n+    it "satisfies `decodeFloat` invariant 2 (zero)" $+      decodeFloat (0 :: BasicDecimal) `shouldBe` (0,0)+    it "satisfies `decodeFloat` invariant 2 (negative zero)" $+      decodeFloat (read "-0" :: BasicDecimal) `shouldBe` (0,0)+    it "satisfies `decodeFloat` invariant 2 (nonzero)" $+      property $ \x ->+      isFinite x && x /= 0 ==> let b      = floatRadix (x :: BasicDecimal)+                                   (m, _) = decodeFloat x+                                   d      = floatDigits x+                                   am     = abs m+                               in b^(d-1) <= am && am < b^d+    it "satisfies relation between `encodeFloat` and `decodeFloat`" $+      property $ \x ->+      not (isNegativeZero x) ==>+      uncurry encodeFloat (decodeFloat x) == (x :: BasicDecimal)+    it "satisfies `exponent` invariant 1 (zero)" $+      exponent (0 :: BasicDecimal) `shouldBe` 0+    it "satisfies `exponent` invariant 1 (nonzero)" $+      property $ \x ->+      isFinite x && x /= 0 ==>+      exponent (x :: BasicDecimal) == snd (decodeFloat x) + floatDigits x+    it "satisfies `exponent` invariant 2" $+      property $ \x ->+      isFinite x ==> let b = floatRadix (x :: BasicDecimal)+                     in x == significand x * fromInteger b ^^ exponent x+    it "satisfies `significand` invariant" $+      property $ \x ->+      isFinite x ==> let s = significand (x :: BasicDecimal)+                         b = floatRadix x+                     in s == 0 ||+                        (s > -1 && s < 1 && abs s >= 1 / fromInteger b)++    it "detects negative zero" $+      isNegativeZero (read "-0" :: BasicDecimal) `shouldBe` True+    it "does not detect normal zero as negative" $+      isNegativeZero (read "+0" :: BasicDecimal) `shouldBe` False+    it "does not detect nonzero numbers as negative zero" $+      property $ \x ->+      x /= 0 ==> isNegativeZero (x :: BasicDecimal) == False++isFinite :: (FinitePrecision p, Rounding r) => Decimal p r -> Bool+isFinite x = not (isNaN x || isInfinite x)

@@ -0,0 +1,618 @@++module Numeric.Decimal.OperationSpec (spec) where++import Test.Hspec++import Numeric.Decimal+import Numeric.Decimal.Arithmetic++import qualified Numeric.Decimal.Operation as Op++spec :: Spec+spec = do+  describe "Special values" $ do+    it "Infinity + 1 = Infinity" $+      op2 Op.add "Infinity" "1" `shouldBe` "Infinity"+    it "NaN + 1 = NaN" $+      op2 Op.add "NaN" "1" `shouldBe` "NaN"+    it "NaN + Infinity = NaN" $+      op2 Op.add "NaN" "Infinity" `shouldBe` "NaN"+    it "1 - Infinity = -Infinity" $+      op2 Op.subtract "1" "Infinity" `shouldBe` "-Infinity"+    it "-1 * Infinity = -Infinity" $+      op2 Op.multiply "-1" "Infinity" `shouldBe` "-Infinity"+    it "-0 - 0 = -0" $+      op2 Op.subtract "-0" "0" `shouldBe` "-0"+    it "-1 * 0 = -0" $+      op2 Op.multiply "-1" "0" `shouldBe` "-0"+    it "1 / 0 = Infinity" $+      op2 Op.divide "1" "0" `shouldBe` "Infinity"+    it "1 / -0 = -Infinity" $+      op2 Op.divide "1" "-0" `shouldBe` "-Infinity"+    it "-1 / 0 = -Infinity" $+      op2 Op.divide "-1" "0" `shouldBe` "-Infinity"++  describe "add" $ do+    it "add('12', '7.00')        ==>  '19.00'" $+      op2 Op.add "12" "7.00" `shouldBe` "19.00"+    it "add('1E+2', '1E+4')      ==>  '1.01E+4'" $+      op2 Op.add "1E+2" "1E+4" `shouldBe` "1.01E+4"++  describe "subtract" $ do+    it "subtract('1.3', '1.07')  ==>  '0.23'" $+      op2 Op.subtract "1.3" "1.07" `shouldBe` "0.23"+    it "subtract('1.3', '1.30')  ==>  '0.00'" $+      op2 Op.subtract "1.3" "1.30" `shouldBe` "0.00"+    it "subtract('1.3', '2.07')  ==>  '-0.77'" $+      op2 Op.subtract "1.3" "2.07" `shouldBe` "-0.77"++  describe "minus" $ do+    it "minus('1.3')   ==>  '-1.3'" $+      op1 Op.minus "1.3" `shouldBe` "-1.3"+    it "minus('-1.3')  ==>  '1.3'" $+      op1 Op.minus "-1.3" `shouldBe` "1.3"++  describe "plus" $ do+    it "plus('1.3')    ==>  '1.3'" $+      op1 Op.plus "1.3" `shouldBe` "1.3"+    it "plus('-1.3')   ==>  '-1.3'" $+      op1 Op.plus "-1.3" `shouldBe` "-1.3"++  describe "multiply" $ do+    it "multiply('1.20', '3')         ==>  '3.60'" $+      op2 Op.multiply "1.20" "3" `shouldBe` "3.60"+    it "multiply('7', '3')            ==>  '21'" $+      op2 Op.multiply "7" "3" `shouldBe` "21"+    it "multiply('0.9', '0.8')        ==> '0.72'" $+      op2 Op.multiply "0.9" "0.8" `shouldBe` "0.72"+    it "multiply('0.9', '-0')         ==> '-0.0'" $+      op2 Op.multiply "0.9" "-0" `shouldBe` "-0.0"+    it "multiply('654321', '654321')  ==>  '4.28135971E+11'" $+      op2 Op.multiply "654321" "654321" `shouldBe` "4.28135971E+11"++  describe "exp" $ do+    it "exp('-Infinity')    ==> '0'" $+      op1 Op.exp "-Infinity" `shouldBe` "0"+    it "exp('-1')           ==> '0.367879441'" $+      op1 Op.exp "-1" `shouldBe` "0.367879441"+    it "exp('0')            ==> '1'" $+      op1 Op.exp "0" `shouldBe` "1"+    it "exp('1')            ==> '2.71828183'" $+      op1 Op.exp "1" `shouldBe` "2.71828183"+    it "exp('0.693147181')  ==> '2.00000000'" $+      op1 Op.exp "0.693147181" `shouldBe` "2.00000000"+    it "exp('+Infinity')    ==> 'Infinity'" $+      op1 Op.exp "+Infinity" `shouldBe` "Infinity"++  describe "fusedMultiplyAdd" $ do+    it ("fused-multiply-add('3', '5', '7')                             ==>  " +++        "'22'") $+      op3 Op.fusedMultiplyAdd "3" "5" "7" `shouldBe` "22"+    it ("fused-multiply-add('3', '-5', '7')                            ==>  " +++        "'-8'") $+      op3 Op.fusedMultiplyAdd "3" "-5" "7" `shouldBe` "-8"+    it ("fused-multiply-add('888565290', '1557.96930', '-86087.7578')  ==>  " +++        "'1.38435736E+12'") $+      op3 Op.fusedMultiplyAdd "888565290" "1557.96930" "-86087.7578" `shouldBe`+      "1.38435736E+12"++  describe "ln" $ do+    it "ln('0')           ==> '-Infinity'" $+      op1 Op.ln "0" `shouldBe` "-Infinity"+    it "ln('1.000')       ==> '0'" $+      op1 Op.ln "1.000" `shouldBe` "0"+    it "ln('2.71828183')  ==> '1.00000000'" $+      op1 Op.ln "2.71828183" `shouldBe` "1.00000000"+    it "ln('10')          ==> '2.30258509'" $+      op1 Op.ln "10" `shouldBe` "2.30258509"+    it "ln('+Infinity')   ==> 'Infinity'" $+      op1 Op.ln "+Infinity" `shouldBe` "Infinity"++  describe "log10" $ do+    it "log10('0')          ==>  '-Infinity'" $+      op1 Op.log10 "0" `shouldBe` "-Infinity"+    it "log10('0.001')      ==>  '-3'" $+      op1 Op.log10 "0.001" `shouldBe` "-3"+    it "log10('1.000')      ==>  '0'" $+      op1 Op.log10 "1.000" `shouldBe` "0"+    it "log10('2')          ==>  '0.301029996'" $+      op1 Op.log10 "2" `shouldBe` "0.301029996"+    it "log10('10')         ==>  '1'" $+      op1 Op.log10 "10" `shouldBe` "1"+    it "log10('70')         ==>  '1.84509804'" $+      op1 Op.log10 "70" `shouldBe` "1.84509804"+    it "log10('+Infinity')  ==>  'Infinity'" $+      op1 Op.log10 "+Infinity" `shouldBe` "Infinity"++  describe "divide" $ do+    it "divide('1', '3'  )      ==>  '0.333333333'" $+      op2 Op.divide "1" "3" `shouldBe` "0.333333333"+    it "divide('2', '3'  )      ==>  '0.666666667'" $+      op2 Op.divide "2" "3" `shouldBe` "0.666666667"+    it "divide('5', '2'  )      ==>  '2.5'" $+      op2 Op.divide "5" "2" `shouldBe` "2.5"+    it "divide('1', '10' )      ==>  '0.1'" $+      op2 Op.divide "1" "10" `shouldBe` "0.1"+    it "divide('12', '12')      ==>  '1'" $+      op2 Op.divide "12" "12" `shouldBe` "1"+    it "divide('8.00', '2')     ==>  '4.00'" $+      op2 Op.divide "8.00" "2" `shouldBe` "4.00"+    it "divide('2.400', '2.0')  ==>  '1.20'" $+      op2 Op.divide "2.400" "2.0" `shouldBe` "1.20"+    it "divide('1000', '100')   ==>  '10'" $+      op2 Op.divide "1000" "100" `shouldBe` "10"+    it "divide('1000', '1')     ==>  '1000'" $+      op2 Op.divide "1000" "1" `shouldBe` "1000"+    it "divide('2.40E+6', '2')  ==>  '1.20E+6'" $+      op2 Op.divide "2.40E+6" "2" `shouldBe` "1.20E+6"++  describe "abs" $ do+    it "abs('2.1')    ==>  '2.1'" $+      op1 Op.abs "2.1" `shouldBe` "2.1"+    it "abs('-100')   ==>  '100'" $+      op1 Op.abs "-100" `shouldBe` "100"+    it "abs('101.5')  ==>  '101.5'" $+      op1 Op.abs "101.5" `shouldBe` "101.5"+    it "abs('-101.5') ==>  '101.5'" $+      op1 Op.abs "-101.5" `shouldBe` "101.5"++  describe "compare" $ do+    it "compare('2.1', '3')     ==>  '-1'" $+      op2 compare' "2.1" "3" `shouldBe` "-1"+    it "compare('2.1', '2.1')   ==>  '0'" $+      op2 compare' "2.1" "2.1" `shouldBe` "0"+    it "compare('2.1', '2.10')  ==>  '0'" $+      op2 compare' "2.1" "2.10" `shouldBe` "0"+    it "compare('3', '2.1')     ==>  '1'" $+      op2 compare' "3" "2.1" `shouldBe` "1"+    it "compare('2.1', '-3')    ==>  '1'" $+      op2 compare' "2.1" "-3" `shouldBe` "1"+    it "compare('-3', '2.1')    ==>  '-1'" $+      op2 compare' "-3" "2.1" `shouldBe` "-1"++  describe "max" $ do+    it "max('3', '2')    ==>  '3'" $+      op2 Op.max "3" "2" `shouldBe` "3"+    it "max('-10', '3')  ==>  '3'" $+      op2 Op.max "-10" "3" `shouldBe` "3"+    it "max('1.0', '1')  ==>  '1'" $+      op2 Op.max "1.0" "1" `shouldBe` "1"+    it "max('7', 'NaN')  ==>  '7'" $+      op2 Op.max "7" "NaN" `shouldBe` "7"++  describe "min" $ do+    it "min('3', '2')    ==>  '2'" $+      op2 Op.min "3" "2" `shouldBe` "2"+    it "min('-10', '3')  ==>  '-10'" $+      op2 Op.min "-10" "3" `shouldBe` "-10"+    it "min('1.0', '1')  ==>  '1.0'" $+      op2 Op.min "1.0" "1" `shouldBe` "1.0"+    it "min('7', 'NaN')  ==>  '7'" $+      op2 Op.min "7" "NaN" `shouldBe` "7"++  describe "power" $ do+    it "power('2', '3')             ==>  '8'" $+      op2 Op.power "2" "3" `shouldBe` "8"+    it "power('-2', '3')            ==>  '-8'" $+      op2 Op.power "-2" "3" `shouldBe` "-8"+    it "power('2', '-3')            ==>  '0.125'" $+      op2 Op.power "2" "-3" `shouldBe` "0.125"+    it "power('1.7', '8')           ==>  '69.7575744'" $+      op2 Op.power "1.7" "8" `shouldBe` "69.7575744"+    it "power('10', '0.301029996')  ==>  '2.00000000'" $+      op2 Op.power "10" "0.301029996" `shouldBe` "2.00000000"+    it "power('Infinity', '-1')     ==>  '0'" $+      op2 Op.power "Infinity" "-1" `shouldBe` "0"+    it "power('Infinity', '0')      ==>  '1'" $+      op2 Op.power "Infinity" "0" `shouldBe` "1"+    it "power('Infinity', '1')      ==>  'Infinity'" $+      op2 Op.power "Infinity" "1" `shouldBe` "Infinity"+    it "power('-Infinity', '-1')    ==>  '-0'" $+      op2 Op.power "-Infinity" "-1" `shouldBe` "-0"+    it "power('-Infinity', '0')     ==>  '1'" $+      op2 Op.power "-Infinity" "0" `shouldBe` "1"+    it "power('-Infinity', '1')     ==>  '-Infinity'" $+      op2 Op.power "-Infinity" "1" `shouldBe` "-Infinity"+    it "power('-Infinity', '2')     ==>  'Infinity'" $+      op2 Op.power "-Infinity" "2" `shouldBe` "Infinity"+    it "power('0', '0')             ==>  'NaN'" $+      op2 Op.power "0" "0" `shouldBe` "NaN"++  describe "quantize" $ do+    it "quantize('2.17', '0.001')        ==>  '2.170'" $+      op2 Op.quantize "2.17" "0.001" `shouldBe` "2.170"+    it "quantize('2.17', '0.01')         ==>  '2.17'" $+      op2 Op.quantize "2.17" "0.01" `shouldBe` "2.17"+    it "quantize('2.17', '0.1')          ==>  '2.2'" $+      op2 Op.quantize "2.17" "0.1" `shouldBe` "2.2"+    it "quantize('2.17', '1e+0')         ==>  '2'" $+      op2 Op.quantize "2.17" "1e+0" `shouldBe` "2"+    it "quantize('2.17', '1e+1')         ==>  '0E+1'" $+      op2 Op.quantize "2.17" "1e+1" `shouldBe` "0E+1"+    it "quantize('-Inf'  'Infinity')     ==>  '-Infinity'" $+      op2 Op.quantize "-Inf" "Infinity" `shouldBe` "-Infinity"+    it "quantize('2',    'Infinity')     ==>  'NaN'" $+      op2 Op.quantize "2" "Infinity" `shouldBe` "NaN"+    it "quantize('-0.1', '1'  )          ==>  '-0'" $+      op2 Op.quantize "-0.1" "1" `shouldBe` "-0"+    it "quantize('-0',   '1e+5')         ==>  '-0E+5'" $+      op2 Op.quantize "-0" "1e+5" `shouldBe` "-0E+5"+    it "quantize('+35236450.6', '1e-2')  ==>  'NaN'" $+      op2 Op.quantize "+35236450.6" "1e-2" `shouldBe` "NaN"+    it "quantize('-35236450.6', '1e-2')  ==>  'NaN'" $+      op2 Op.quantize "-35236450.6" "1e-2" `shouldBe` "NaN"+    it "quantize('217',  '1e-1')         ==>  '217.0'" $+      op2 Op.quantize "217" "1e-1" `shouldBe` "217.0"+    it "quantize('217',  '1e+0')         ==>  '217'" $+      op2 Op.quantize "217" "1e+0" `shouldBe` "217"+    it "quantize('217',  '1e+1')         ==>  '2.2E+2'" $+      op2 Op.quantize "217" "1e+1" `shouldBe` "2.2E+2"+    it "quantize('217',  '1e+2')         ==>  '2E+2'" $+      op2 Op.quantize "217" "1e+2" `shouldBe` "2E+2"++  describe "reduce" $ do+    it "reduce('2.1')     ==>  '2.1'" $+      op1 Op.reduce "2.1" `shouldBe` "2.1"+    it "reduce('-2.0')    ==>  '-2'" $+      op1 Op.reduce "-2.0" `shouldBe` "-2"+    it "reduce('1.200')   ==>  '1.2'" $+      op1 Op.reduce "1.200" `shouldBe` "1.2"+    it "reduce('-120')    ==>  '-1.2E+2'" $+      op1 Op.reduce "-120" `shouldBe` "-1.2E+2"+    it "reduce('120.00')  ==>  '1.2E+2'" $+      op1 Op.reduce "120.00" `shouldBe` "1.2E+2"+    it "reduce('0.00')    ==>  '0'" $+      op1 Op.reduce "0.00" `shouldBe` "0"++  describe "roundToIntegralExact" $ do+    it "round-to-integral-exact('2.1')      ==>  '2'" $+      op1 Op.roundToIntegralExact "2.1" `shouldBe` "2"+    it "round-to-integral-exact('100')      ==>  '100'" $+      op1 Op.roundToIntegralExact "100" `shouldBe` "100"+    it "round-to-integral-exact('100.0')    ==>  '100'" $+      op1 Op.roundToIntegralExact "100.0" `shouldBe` "100"+    it "round-to-integral-exact('101.5')    ==>  '102'" $+      op1 Op.roundToIntegralExact "101.5" `shouldBe` "102"+    it "round-to-integral-exact('-101.5')   ==>  '-102'" $+      op1 Op.roundToIntegralExact "-101.5" `shouldBe` "-102"+    it "round-to-integral-exact('10E+5')    ==>  '1.0E+6'" $+      op1 Op.roundToIntegralExact "10E+5" `shouldBe` "1.0E+6"+    it "round-to-integral-exact('7.89E+77') ==>  '7.89E+77'" $+      op1 Op.roundToIntegralExact "7.89E+77" `shouldBe` "7.89E+77"+    it "round-to-integral-exact('-Inf')     ==>  '-Infinity'" $+      op1 Op.roundToIntegralExact "-Inf" `shouldBe` "-Infinity"++  describe "squareRoot" $ do+    it "square-root('0')     ==> '0'" $+      op1 Op.squareRoot "0" `shouldBe` "0"+    it "square-root('-0')    ==> '-0'" $+      op1 Op.squareRoot "-0" `shouldBe` "-0"+    -- The following example is a corrected version of that found in the+    -- specification; confirmed with Mike Cowlishaw on 2016-08-02.+    it "square-root('0.39')  ==> '0.624499800'" $+      op1 Op.squareRoot "0.39" `shouldBe` "0.624499800"+    it "square-root('100')   ==> '10'" $+      op1 Op.squareRoot "100" `shouldBe` "10"+    it "square-root('1')     ==> '1'" $+      op1 Op.squareRoot "1" `shouldBe` "1"+    it "square-root('1.0')   ==> '1.0'" $+      op1 Op.squareRoot "1.0" `shouldBe` "1.0"+    it "square-root('1.00')  ==> '1.0'" $+      op1 Op.squareRoot "1.00" `shouldBe` "1.0"+    it "square-root('7')     ==> '2.64575131'" $+      op1 Op.squareRoot "7" `shouldBe` "2.64575131"+    it "square-root('10')    ==> '3.16227766'" $+      op1 Op.squareRoot "10" `shouldBe` "3.16227766"++  describe "and" $ do+    it "and('0', '0')        ==>  '0'" $+      op2 Op.and "0" "0" `shouldBe` "0"+    it "and('0', '1')        ==>  '0'" $+      op2 Op.and "0" "1" `shouldBe` "0"+    it "and('1', '0')        ==>  '0'" $+      op2 Op.and "1" "0" `shouldBe` "0"+    it "and('1', '1')        ==>  '1'" $+      op2 Op.and "1" "1" `shouldBe` "1"+    it "and('1100', '1010')  ==>  '1000'" $+      op2 Op.and "1100" "1010" `shouldBe` "1000"+    it "and('1111', '10')    ==>  '10'" $+      op2 Op.and "1111" "10" `shouldBe` "10"++  describe "or" $ do+    it "or('0', '0')        ==>  '0'" $+      op2 Op.or "0" "0" `shouldBe` "0"+    it "or('0', '1')        ==>  '1'" $+      op2 Op.or "0" "1" `shouldBe` "1"+    it "or('1', '0')        ==>  '1'" $+      op2 Op.or "1" "0" `shouldBe` "1"+    it "or('1', '1')        ==>  '1'" $+      op2 Op.or "1" "1" `shouldBe` "1"+    it "or('1100', '1010')  ==>  '1110'" $+      op2 Op.or "1100" "1010" `shouldBe` "1110"+    it "or('1110', '10')    ==>  '1110'" $+      op2 Op.or "1110" "10" `shouldBe` "1110"++  describe "xor" $ do+    it "xor('0', '0')        ==>  '0'" $+      op2 Op.xor "0" "0" `shouldBe` "0"+    it "xor('0', '1')        ==>  '1'" $+      op2 Op.xor "0" "1" `shouldBe` "1"+    it "xor('1', '0')        ==>  '1'" $+      op2 Op.xor "1" "0" `shouldBe` "1"+    it "xor('1', '1')        ==>  '0'" $+      op2 Op.xor "1" "1" `shouldBe` "0"+    it "xor('1100', '1010')  ==>  '110'" $+      op2 Op.xor "1100" "1010" `shouldBe` "110"+    it "xor('1111', '10')    ==>  '1101'" $+      op2 Op.xor "1111" "10" `shouldBe` "1101"++  describe "invert" $ do+    it "invert('0')          ==>  '111111111'" $+      op1 Op.invert "0" `shouldBe` "111111111"+    it "invert('1')          ==>  '111111110'" $+      op1 Op.invert "1" `shouldBe` "111111110"+    it "invert('111111111')  ==>  '0'" $+      op1 Op.invert "111111111" `shouldBe` "0"+    it "invert('101010101')  ==>  '10101010'" $+      op1 Op.invert "101010101" `shouldBe` "10101010"++  describe "canonical" $+    it "canonical('2.50')  ==>  '2.50'" $+      op1 Op.canonical "2.50" `shouldBe` "2.50"++  describe "class'" $ do+    it "class('Infinity')   ==>  \"+Infinity\"" $+      op1 Op.class' "Infinity" `shouldBe` "+Infinity"+    it "class('1E-10')      ==>  \"+Normal\"" $+      op1 Op.class' "1E-10" `shouldBe` "+Normal"+    it "class('2.50')       ==>  \"+Normal\"" $+      op1 Op.class' "2.50" `shouldBe` "+Normal"+    it "class('0.1E-999')   ==>  \"+Subnormal\"" $+      op1 Op.class' "0.1E-999" `shouldBe` "+Subnormal"+    it "class('0')          ==>  \"+Zero\"" $+      op1 Op.class' "0" `shouldBe` "+Zero"+    it "class('-0')         ==>  \"-Zero\"" $+      op1 Op.class' "-0" `shouldBe` "-Zero"+    it "class('-0.1E-999')  ==>  \"-Subnormal\"" $+      op1 Op.class' "-0.1E-999" `shouldBe` "-Subnormal"+    it "class('-1E-10')     ==>  \"-Normal\"" $+      op1 Op.class' "-1E-10" `shouldBe` "-Normal"+    it "class('-2.50')      ==>  \"-Normal\"" $+      op1 Op.class' "-2.50" `shouldBe` "-Normal"+    it "class('-Infinity')  ==>  \"-Infinity\"" $+      op1 Op.class' "-Infinity" `shouldBe` "-Infinity"+    it "class('NaN')        ==>  \"NaN\"" $+      op1 Op.class' "NaN" `shouldBe` "NaN"+    it "class('-NaN')       ==>  \"NaN\"" $+      op1 Op.class' "-NaN" `shouldBe` "NaN"+    it "class('sNaN')       ==>  \"sNaN\"" $+      op1 Op.class' "sNaN" `shouldBe` "sNaN"++  describe "compareTotal" $ do+    it "compare-total('12.73', '127.9')   ==>  '-1'" $+      op2 compareTotal' "12.73" "127.9" `shouldBe` "-1"+    it "compare-total('-127',  '12')      ==>  '-1'" $+      op2 compareTotal' "-127" "12" `shouldBe` "-1"+    it "compare-total('12.30', '12.3')    ==>  '-1'" $+      op2 compareTotal' "12.30" "12.3" `shouldBe` "-1"+    it "compare-total('12.30', '12.30')   ==>  '0'" $+      op2 compareTotal' "12.30" "12.30" `shouldBe` "0"+    it "compare-total('12.3',  '12.300')  ==>  '1'" $+      op2 compareTotal' "12.3" "12.300" `shouldBe` "1"+    it "compare-total('12.3',  'NaN')     ==>  '-1'" $+      op2 compareTotal' "12.3" "NaN" `shouldBe` "-1"++  describe "copy" $ do+    it "copy('2.1')    ==>  '2.1'" $+      op1 Op.copy "2.1" `shouldBe` "2.1"+    it "copy('-1.00')  ==>  '-1.00'" $+      op1 Op.copy "-1.00" `shouldBe` "-1.00"++  describe "copyAbs" $ do+    it "copy-abs('2.1')   ==>  '2.1'" $+      op1 Op.copyAbs "2.1" `shouldBe` "2.1"+    it "copy-abs('-100')  ==>  '100'" $+      op1 Op.copyAbs "-100" `shouldBe` "100"++  describe "copyNegate" $ do+    it "copy-negate('101.5')   ==>  '-101.5'" $+      op1 Op.copyNegate "101.5" `shouldBe` "-101.5"+    it "copy-negate('-101.5')  ==>  '101.5'" $+      op1 Op.copyNegate "-101.5" `shouldBe` "101.5"++  describe "copySign" $ do+    it "copy-sign( '1.50',  '7.33')  ==>  '1.50'" $+      op2 Op.copySign  "1.50"  "7.33" `shouldBe` "1.50"+    it "copy-sign('-1.50',  '7.33')  ==>  '1.50'" $+      op2 Op.copySign "-1.50"  "7.33" `shouldBe` "1.50"+    it "copy-sign( '1.50', '-7.33')  ==>  '-1.50'" $+      op2 Op.copySign  "1.50" "-7.33" `shouldBe` "-1.50"+    it "copy-sign('-1.50', '-7.33')  ==>  '-1.50'" $+      op2 Op.copySign "-1.50" "-7.33" `shouldBe` "-1.50"++  describe "isCanonical" $+    it "is-canonical('2.50')  ==>  '1'" $+      pred1 Op.isCanonical "2.50" `shouldBe` "1"++  describe "isFinite" $ do+    it "is-finite('2.50')  ==>  '1'" $+      pred1 Op.isFinite "2.50" `shouldBe` "1"+    it "is-finite('-0.3')  ==>  '1'" $+      pred1 Op.isFinite "-0.3" `shouldBe` "1"+    it "is-finite('0')     ==>  '1'" $+      pred1 Op.isFinite "0" `shouldBe` "1"+    it "is-finite('Inf')   ==>  '0'" $+      pred1 Op.isFinite "Inf" `shouldBe` "0"+    it "is-finite('NaN')   ==>  '0'" $+      pred1 Op.isFinite "NaN" `shouldBe` "0"++  describe "isInfinite" $ do+    it "is-infinite('2.50')  ==>  '0'" $+      pred1 Op.isInfinite "2.50" `shouldBe` "0"+    it "is-infinite('-Inf')  ==>  '1'" $+      pred1 Op.isInfinite "-Inf" `shouldBe` "1"+    it "is-infinite('NaN')   ==>  '0'" $+      pred1 Op.isInfinite "NaN" `shouldBe` "0"++  describe "isNaN" $ do+    it "is-NaN('2.50')   ==>  '0'" $+      pred1 Op.isNaN "2.50" `shouldBe` "0"+    it "is-NaN('NaN')    ==>  '1'" $+      pred1 Op.isNaN "NaN" `shouldBe` "1"+    it "is-NaN('-sNaN')  ==>  '1'" $+      pred1 Op.isNaN "-sNaN" `shouldBe` "1"++  describe "isNormal" $ do+    it "is-normal('2.50')      ==>  '1'" $+      pred1 Op.isNormal "2.50" `shouldBe` "1"+    it "is-normal('0.1E-999')  ==>  '0'" $+      pred1 Op.isNormal "0.1E-999" `shouldBe` "0"+    it "is-normal('0.00')      ==>  '0'" $+      pred1 Op.isNormal "0.00" `shouldBe` "0"+    it "is-normal('-Inf')      ==>  '0'" $+      pred1 Op.isNormal "-Inf" `shouldBe` "0"+    it "is-normal('NaN')       ==>  '0'" $+      pred1 Op.isNormal "NaN" `shouldBe` "0"++  describe "isQNaN" $ do+    it "is-qNaN('2.50')  ==>  '0'" $+      pred1 Op.isQNaN "2.50" `shouldBe` "0"+    it "is-qNaN('NaN')   ==>  '1'" $+      pred1 Op.isQNaN "NaN" `shouldBe` "1"+    it "is-qNaN('sNaN')  ==>  '0'" $+      pred1 Op.isQNaN "sNaN" `shouldBe` "0"++  describe "isSigned" $ do+    it "is-signed('2.50')  ==>  '0'" $+      pred1 Op.isSigned "2.50" `shouldBe` "0"+    it "is-signed('-12')   ==>  '1'" $+      pred1 Op.isSigned "-12" `shouldBe` "1"+    it "is-signed('-0')    ==>  '1'" $+      pred1 Op.isSigned "-0" `shouldBe` "1"++  describe "isSNaN" $ do+    it "is-sNaN('2.50')  ==>  '0'" $+      pred1 Op.isSNaN "2.50" `shouldBe` "0"+    it "is-sNaN('NaN')   ==>  '0'" $+      pred1 Op.isSNaN "NaN" `shouldBe` "0"+    it "is-sNaN('sNaN')  ==>  '1'" $+      pred1 Op.isSNaN "sNaN" `shouldBe` "1"++  describe "isSubnormal" $ do+    it "is-subnormal('2.50')      ==>  '0'" $+      pred1 Op.isSubnormal "2.50" `shouldBe` "0"+    it "is-subnormal('0.1E-999')  ==>  '1'" $+      pred1 Op.isSubnormal "0.1E-999" `shouldBe` "1"+    it "is-subnormal('0.00')      ==>  '0'" $+      pred1 Op.isSubnormal "0.00" `shouldBe` "0"+    it "is-subnormal('-Inf')      ==>  '0'" $+      pred1 Op.isSubnormal "-Inf" `shouldBe` "0"+    it "is-subnormal('NaN')       ==>  '0'" $+      pred1 Op.isSubnormal "NaN" `shouldBe` "0"++  describe "isZero" $ do+    it "is-zero('0')      ==>  '1'" $+      pred1 Op.isZero "0" `shouldBe` "1"+    it "is-zero('2.50')   ==>  '0'" $+      pred1 Op.isZero "2.50" `shouldBe` "0"+    it "is-zero('-0E+2')  ==>  '1'" $+      pred1 Op.isZero "-0E+2" `shouldBe` "1"++  describe "logb" $ do+    it "logb('250')   ==>  '2'" $+      op1 Op.logb "250" `shouldBe` "2"+    it "logb('2.50')  ==>  '0'" $+      op1 Op.logb "2.50" `shouldBe` "0"+    it "logb('0.03')  ==>  '-2'" $+      op1 Op.logb "0.03" `shouldBe` "-2"+    it "logb('0')     ==>  '-Infinity'" $+      op1 Op.logb "0" `shouldBe` "-Infinity"++  describe "scaleb" $ do+    it "scaleb('7.50', '-2')  ==>  '0.0750'" $+      op2 Op.scaleb "7.50" "-2" `shouldBe` "0.0750"+    it "scaleb('7.50', '0')   ==>  '7.50'" $+      op2 Op.scaleb "7.50" "0" `shouldBe` "7.50"+    it "scaleb('7.50', '3')   ==>  '7.50E+3'" $+      op2 Op.scaleb "7.50" "3" `shouldBe` "7.50E+3"++  describe "radix" $+    it "radix()  ==>  '10'" $+      op0 Op.radix `shouldBe` "10"++  describe "sameQuantum" $ do+    it "samequantum('2.17', '0.001')  ==>  '0'" $+      pred2 Op.sameQuantum "2.17" "0.001" `shouldBe` "0"+    it "samequantum('2.17', '0.01')   ==>  '1'" $+      pred2 Op.sameQuantum "2.17" "0.01" `shouldBe` "1"+    it "samequantum('2.17', '0.1')    ==>  '0'" $+      pred2 Op.sameQuantum "2.17" "0.1" `shouldBe` "0"+    it "samequantum('2.17', '1')      ==>  '0'" $+      pred2 Op.sameQuantum "2.17" "1" `shouldBe` "0"+    it "samequantum('Inf', '-Inf')    ==>  '1'" $+      pred2 Op.sameQuantum "Inf" "-Inf" `shouldBe` "1"+    it "samequantum('NaN', 'NaN')     ==>  '1'" $+      pred2 Op.sameQuantum "NaN" "NaN" `shouldBe` "1"++  describe "shift" $ do+    it "shift('34', '8')          ==>  '400000000'" $+      op2 Op.shift "34" "8" `shouldBe` "400000000"+    it "shift('12', '9')          ==>  '0'" $+      op2 Op.shift "12" "9" `shouldBe` "0"+    it "shift('123456789', '-2')  ==>  '1234567'" $+      op2 Op.shift "123456789" "-2" `shouldBe` "1234567"+    it "shift('123456789', '0')   ==>  '123456789'" $+      op2 Op.shift "123456789" "0" `shouldBe` "123456789"+    it "shift('123456789', '+2')  ==>  '345678900'" $+      op2 Op.shift "123456789" "+2" `shouldBe` "345678900"++  describe "rotate" $ do+    it "rotate('34', '8')          ==>  '400000003'" $+      op2 Op.rotate "34" "8" `shouldBe` "400000003"+    it "rotate('12', '9')          ==>  '12'" $+      op2 Op.rotate "12" "9" `shouldBe` "12"+    it "rotate('123456789', '-2')  ==>  '891234567'" $+      op2 Op.rotate "123456789" "-2" `shouldBe` "891234567"+    it "rotate('123456789', '0')   ==>  '123456789'" $+      op2 Op.rotate "123456789" "0" `shouldBe` "123456789"+    it "rotate('123456789', '+2')  ==>  '345678912'" $+      op2 Op.rotate "123456789" "+2" `shouldBe` "345678912"++exceptionError :: Exception p r -> a+exceptionError = error . show . exceptionSignal++type BasicArith = Arith P9 RoundHalfUp++pred1 :: (BasicDecimal -> BasicArith Bool) -> String -> String+pred1 op x = either exceptionError (show . fromBool) $+  evalArith arith newContext+  where arith = op (read x)++pred2 :: (BasicDecimal -> BasicDecimal -> BasicArith Bool) -> String -> String+      -> String+pred2 op x y = either exceptionError (show . fromBool) $+  evalArith arith newContext+  where arith = op (read x) (read y)++op0 :: Show a => BasicArith a -> String+op0 op = either exceptionError show $ evalArith op newContext++op1 :: Show a => (BasicDecimal -> BasicArith a) -> String -> String+op1 op x = either exceptionError show $ evalArith arith newContext+  where arith = op (read x)++op2 :: Show a => (BasicDecimal -> BasicDecimal -> BasicArith a) -> String+    -> String -> String+op2 op x y = either exceptionError show $ evalArith arith newContext+  where arith = read x `op` read y++op3 :: Show a => (BasicDecimal -> BasicDecimal -> BasicDecimal -> BasicArith a)+    -> String -> String -> String -> String+op3 op x y z = either exceptionError show $ evalArith arith newContext+  where arith = op (read x) (read y) (read z)++compare' :: Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+compare' x y = either id fromOrdering <$> Op.compare x y++compareTotal' :: Decimal a b -> Decimal c d -> Arith p r (Decimal p r)+compareTotal' x y = fromOrdering <$> Op.compareTotal x y

@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}