diff --git a/scientific.cabal b/scientific.cabal --- a/scientific.cabal +++ b/scientific.cabal @@ -1,5 +1,5 @@ name: scientific -version: 0.3.0.2 +version: 0.3.1.0 synopsis: Numbers represented using scientific notation description: @Data.Scientific@ provides a space efficient and arbitrary precision @@ -18,8 +18,7 @@ . * A 'Scientific' is more efficient to construct. Rational numbers need to be constructed using '%' which has to compute the 'gcd' of the 'numerator' and - 'denominator'. Scientific numbers only need to be normalized, i.e. @10000000@ - to @1e7@. + 'denominator'. . * 'Scientific' is safe against numbers with huge exponents. For example: @1e1000000000 :: 'Rational'@ will fill up all space and crash your @@ -50,15 +49,17 @@ exposed-modules: Data.Scientific Data.Text.Lazy.Builder.Scientific Data.ByteString.Builder.Scientific + other-modules: Math.NumberTheory.Logarithms other-extensions: DeriveDataTypeable, BangPatterns ghc-options: -Wall - build-depends: base >= 4.3 && < 4.8 - , deepseq >= 1.3 && < 1.4 - , text >= 0.8 && < 1.3 - , bytestring >= 0.10 && < 0.11 - , hashable >= 1.1.2 && < 1.3 - , arithmoi >= 0.4.1 && < 0.5 - , array >= 0.1 && < 0.6 + build-depends: base >= 4.3 && < 4.8 + , ghc-prim + , integer-gmp + , deepseq >= 1.3 && < 1.4 + , text >= 0.8 && < 1.3 + , bytestring >= 0.10 && < 0.11 + , hashable >= 1.1.2 && < 1.3 + , array >= 0.1 && < 0.6 hs-source-dirs: src default-language: Haskell2010 @@ -71,7 +72,8 @@ build-depends: scientific , base >= 4.3 && < 4.8 - , tasty >= 0.3.1 && < 0.9 + , tasty >= 0.5 && < 0.9 + , tasty-ant-xml >= 1.0 && < 1.1 , tasty-smallcheck >= 0.2 && < 0.9 , tasty-quickcheck >= 0.8 && < 0.9 , smallcheck >= 1.0 && < 1.2 @@ -81,10 +83,16 @@ benchmark bench-scientific type: exitcode-stdio-1.0 - hs-source-dirs: bench + hs-source-dirs: bench src main-is: bench.hs default-language: Haskell2010 ghc-options: -O2 - build-depends: scientific - , base >= 4.3 && < 4.8 + build-depends: base >= 4.3 && < 4.8 , criterion >= 0.5 && < 0.12 + , ghc-prim + , integer-gmp + , deepseq >= 1.3 && < 1.4 + , text >= 0.8 && < 1.3 + , bytestring >= 0.10 && < 0.11 + , hashable >= 1.1.2 && < 1.3 + , array >= 0.1 && < 0.6 diff --git a/src/Data/Scientific.hs b/src/Data/Scientific.hs --- a/src/Data/Scientific.hs +++ b/src/Data/Scientific.hs @@ -22,8 +22,7 @@ -- -- * A 'Scientific' is more efficient to construct. Rational numbers need to be -- constructed using '%' which has to compute the 'gcd' of the 'numerator' and --- 'denominator'. Scientific numbers only need to be normalized, i.e. @10000000@ --- to @1e7@. +-- 'denominator'. -- -- * 'Scientific' is safe against numbers with huge exponents. For example: -- @1e1000000000 :: 'Rational'@ will fill up all space and crash your @@ -36,6 +35,11 @@ -- to @'Integral's@ (like: 'Int') or @'RealFloat's@ (like: 'Double' or 'Float') -- will always be bounded by the target type. -- +-- Note that, in order to do fast magnitude computations (@10^e@), this module +-- computes the first 1100 powers of 10 and stores them in a top-level +-- array. This means that the first magnitude computation (used in 'realToFrac' +-- among others) is slower but subsequent computations should be O(1). +-- -- This module is designed to be imported qualified: -- -- @import Data.Scientific as Scientific@ @@ -58,6 +62,9 @@ , FPFormat(..) , toDecimalDigits + + -- * Normalization + , normalize ) where @@ -99,6 +106,14 @@ data Scientific = Scientific { coefficient :: !Integer -- ^ The coefficient of a scientific number. + -- + -- Note that this number is not necessarily normalized, i.e. + -- it could contain trailing zeros. + -- + -- Scientific numbers are automatically normalized when pretty printed or + -- in 'toDecimalDigits'. + -- + -- Use 'normalize' to do manual normalization. , base10Exponent :: {-# UNPACK #-} !Int -- ^ The base-10 exponent of a scientific number. @@ -106,20 +121,8 @@ -- | @scientific c e@ constructs a scientific number which corresponds -- to the 'Fractional' number: @'fromInteger' c * 10 '^^' e@. --- --- Note that this function performs normalization, i.e. it divides out powers of --- 10 from @c@ and adds them to @e@. scientific :: Integer -> Int -> Scientific -scientific c !e - | c > 0 = normalize c e - | c == 0 = Scientific 0 0 - | otherwise = -(normalize (-c) e) -{-# INLINE scientific #-} - -normalize :: Integer -> Int -> Scientific -normalize c !e = case quotRem c 10 of - (q, 0) -> normalize q (e+1) - _ -> Scientific c e +scientific = Scientific ---------------------------------------------------------------------- @@ -318,7 +321,7 @@ -- space usage is bounded by the target type. -- -- For large negative exponents we check if the exponent is smaller --- than some limit (currently -20). In that case we know that the +-- than some limit (currently -1100). In that case we know that the -- scientific number is really small (unless the coefficient has many -- digits) so we can immediately return -1 for negative scientific -- numbers or 0 for positive numbers. @@ -331,11 +334,11 @@ -- (log10 c) if the exponent is not below the limit. dangerouslySmall :: Integer -> Int -> Bool dangerouslySmall c e = e < (-limit) && e < (-integerLog10' (abs c)) - 1 - where - limit :: Int - limit = 20 {-# INLINE dangerouslySmall #-} +limit :: Int +limit = maxExpt + positivize :: (Ord a, Num a, Num b) => (a -> b) -> (a -> b) positivize f x | x < 0 = -(f (-x)) | otherwise = f x @@ -407,27 +410,32 @@ -- scientific numbers coming from an untrusted source. toRealFloat :: forall a. (RealFloat a) => Scientific -> a toRealFloat s@(Scientific c e) - | e > hiLimit = sign (1/0) -- Infinity - | e < loLimit && e + d < loLimit = sign 0 - | otherwise = realToFrac s + | e > limit && e > hiLimit = sign (1/0) -- Infinity + | e < -limit && e < loLimit && e + d < loLimit = sign 0 + | otherwise = realToFrac s where - hiLimit = ceiling (fromIntegral hi * log10Radix) - loLimit = floor (fromIntegral lo * log10Radix) - - ceiling (fromIntegral digits * log10Radix) - - log10Radix :: Double - log10Radix = logBase 10 $ fromInteger radix - - radix = floatRadix (undefined :: a) - digits = floatDigits (undefined :: a) - (lo, hi) = floatRange (undefined :: a) + (loLimit, hiLimit) = exponentLimits (undefined :: a) d = integerLog10' (abs c) sign x | c < 0 = -x | otherwise = x +exponentLimits :: forall a. (RealFloat a) => a -> (Int, Int) +exponentLimits _ = (loLimit, hiLimit) + where + loLimit = floor (fromIntegral lo * log10Radix) - + ceiling (fromIntegral digits * log10Radix) + hiLimit = ceiling (fromIntegral hi * log10Radix) + log10Radix :: Double + log10Radix = logBase 10 $ fromInteger radix + + radix = floatRadix (undefined :: a) + digits = floatDigits (undefined :: a) + (lo, hi) = floatRange (undefined :: a) + + ---------------------------------------------------------------------- -- Parsing ---------------------------------------------------------------------- @@ -589,7 +597,7 @@ ---------------------------------------------------------------------- --- | Similar to 'floatToDigits', @toDecimalDigits@ takes a +-- | Similar to 'Numeric.floatToDigits', @toDecimalDigits@ takes a -- non-negative 'Scientific' number, and returns a list of digits and -- a base-10 exponent. In particular, if @x>=0@, and -- @@ -602,10 +610,17 @@ -- (2) @x = 0.d1d2...dn * (10^^e)@ -- -- (3) @0 <= di <= 9@ +-- +-- (4) @null $ takeWhile (==0) $ reverse [d1,d2,...,dn]@ +-- +-- The last property means that the coefficient will be normalized, i.e. doesn't +-- contain trailing zeros. toDecimalDigits :: Scientific -> ([Int], Int) -toDecimalDigits (Scientific 0 _) = ([0], 0) -toDecimalDigits (Scientific c e) = (is, n + e) +toDecimalDigits (Scientific 0 _) = ([0], 0) +toDecimalDigits (Scientific c' e') = (is, n + e) where + Scientific c e = normalizePositive c' e' + (is, n) = reverseAndLength $ digits c digits :: Integer -> [Int] @@ -619,3 +634,24 @@ where rev [] a !m = (a, m) rev (x:xs) a !m = rev xs (x:a) (m+1) + + +---------------------------------------------------------------------- +-- Normalization +---------------------------------------------------------------------- + +-- | Normalize a scientific number by dividing out powers of 10 from the +-- 'coefficient' and incrementing the 'base10Exponent' each time. +-- +-- You should rarely have a need for this function since scientific numbers are +-- automatically normalized when pretty-printed and in 'toDecimalDigits'. +normalize :: Scientific -> Scientific +normalize (Scientific c e) + | c < 0 = -(normalizePositive (-c) e) + | c > 0 = normalizePositive c e + | otherwise {- c == 0 -} = Scientific 0 0 + +normalizePositive :: Integer -> Int -> Scientific +normalizePositive c !e = case quotRem c 10 of + (q, 0) -> normalizePositive q (e+1) + _ -> Scientific c e diff --git a/src/Math/NumberTheory/Logarithms.hs b/src/Math/NumberTheory/Logarithms.hs new file mode 100644 --- /dev/null +++ b/src/Math/NumberTheory/Logarithms.hs @@ -0,0 +1,60 @@ +{-# LANGUAGE CPP, MagicHash #-} + +-- | Integer logarithm, copied from @arithmoi@ +module Math.NumberTheory.Logarithms ( integerLog10' ) where + +import GHC.Base +#if __GLASGOW_HASKELL__ < 705 +import GHC.Word +#endif +import GHC.Integer.Logarithms + +-- | Only defined for positive inputs! +integerLog10' :: Integer -> Int +integerLog10' n + | n < 10 = 0 + | n < 100 = 1 + | otherwise = ex + integerLog10' (n `quot` integerPower 10 ex) + where + ln = I# (integerLog2# n) + -- u/v is a good approximation of log 2/log 10 + u = 1936274 + v = 6432163 + -- so ex is a good approximation to integerLogBase 10 n + ex = fromInteger ((u * fromIntegral ln) `quot` v) + +-- | Power of an 'Integer' by the left-to-right repeated squaring algorithm. +-- This needs two multiplications in each step while the right-to-left +-- algorithm needs only one multiplication for 0-bits, but here the +-- two factors always have approximately the same size, which on average +-- gains a bit when the result is large. +-- +-- For small results, it is unlikely to be any faster than '(^)', quite +-- possibly slower (though the difference shouldn't be large), and for +-- exponents with few bits set, the same holds. But for exponents with +-- many bits set, the speedup can be significant. +-- +-- /Warning:/ No check for the negativity of the exponent is performed, +-- a negative exponent is interpreted as a large positive exponent. +integerPower :: Integer -> Int -> Integer +integerPower b (I# e#) = power b (int2Word# e#) + +power :: Integer -> Word# -> Integer +power b w# + | isTrue# (w# `eqWord#` 0##) = 1 + | isTrue# (w# `eqWord#` 1##) = b + | otherwise = go (wordLog2# w# -# 1#) b (b*b) + where + go 0# l h = if isTrue# ((w# `and#` 1##) `eqWord#` 0##) then l*l else (l*h) + go i# l h + | w# `hasBit#` i# = go (i# -# 1#) (l*h) (h*h) + | otherwise = go (i# -# 1#) (l*l) (l*h) + +-- | A raw version of testBit for 'Word#'. +hasBit# :: Word# -> Int# -> Bool +hasBit# w# i# = isTrue# (((w# `uncheckedShiftRL#` i#) `and#` 1##) `neWord#` 0##) + +#if __GLASGOW_HASKELL__ < 707 +isTrue# :: Bool -> Bool +isTrue# = id +#endif diff --git a/test/test.hs b/test/test.hs --- a/test/test.hs +++ b/test/test.hs @@ -13,6 +13,7 @@ import Control.Monad import Data.Scientific as Scientific import Test.Tasty +import Test.Tasty.Runners.AntXML import qualified Test.SmallCheck as SC import qualified Test.SmallCheck.Series as SC import qualified Test.Tasty.SmallCheck as SC (testProperty) @@ -31,10 +32,12 @@ #endif main :: IO () -main = defaultMain $ testGroup "scientific" +main = testMain $ testGroup "scientific" [ smallQuick "normalization" - (\s -> s /= 0 SC.==> abs (Scientific.coefficient s) `mod` 10 /= 0) - (\s -> s /= 0 QC.==> abs (Scientific.coefficient s) `mod` 10 /= 0) + (SC.over normalizedScientificSeries $ \s -> + s /= 0 SC.==> abs (Scientific.coefficient s) `mod` 10 /= 0) + (QC.forAll normalizedScientificGen $ \s -> + s /= 0 QC.==> abs (Scientific.coefficient s) `mod` 10 /= 0) , testGroup "Formatting" [ testProperty "read . show == id" $ \s -> read (show s) === s @@ -131,6 +134,9 @@ ] ] +testMain :: TestTree -> IO () +testMain = defaultMainWithIngredients (antXMLRunner:defaultIngredients) + conversionsProperties :: forall realFloat. ( RealFloat realFloat , QC.Arbitrary realFloat @@ -183,15 +189,19 @@ toDecimalDigits_laws x = let (ds, e) = Scientific.toDecimalDigits x - rule1 = length ds >= 1 + rule1 = n >= 1 + n = length ds rule2 = toRational x == coeff * 10 ^^ e coeff = foldr (\di a -> a / 10 + fromIntegral di) 0 (0:ds) rule3 = all (\di -> 0 <= di && di <= 9) ds - in rule1 && rule2 && rule3 + rule4 | n == 1 = True + | otherwise = null $ takeWhile (==0) $ reverse ds + in rule1 && rule2 && rule3 && rule4 + properFraction_laws :: Scientific -> Bool properFraction_laws x = fromInteger n + f === x && (positive n == posX || n == 0) && @@ -238,7 +248,10 @@ nonNegativeScientificSeries :: (Monad m) => SC.Series m Scientific nonNegativeScientificSeries = liftM SC.getNonNegative SC.series +normalizedScientificSeries :: (Monad m) => SC.Series m Scientific +normalizedScientificSeries = liftM Scientific.normalize SC.series + ---------------------------------------------------------------------- -- QuickCheck instances ---------------------------------------------------------------------- @@ -253,6 +266,9 @@ nonNegativeScientificGen = scientific <$> (QC.getNonNegative <$> QC.arbitrary) <*> intGen + +normalizedScientificGen :: QC.Gen Scientific +normalizedScientificGen = Scientific.normalize <$> QC.arbitrary intGen :: QC.Gen Int #if MIN_VERSION_QuickCheck(2,7,0)