diff --git a/scientific.cabal b/scientific.cabal --- a/scientific.cabal +++ b/scientific.cabal @@ -1,18 +1,37 @@ name: scientific -version: 0.2.0.2 -synopsis: Arbitrary-precision floating-point numbers represented using scientific notation -description: A @Scientific@ number is an arbitrary-precision floating-point number - represented using scientific notation. - . - A scientific number with 'coefficient' @c@ and - 'base10Exponent' @e@ corresponds to the - 'Fractional' number: @'fromInteger' c * 10 '^^' e@ - . - Its primary use-case is to serve as the target of - parsing floating point numbers. Since the textual - representation of floating point numbers use - scientific notation they can be efficiently - parsed to a @Scientific@ number. +version: 0.3.0.0 +synopsis: Numbers represented using scientific notation +description: + @Data.Scientific@ provides a space efficient and arbitrary precision + scientific number type. + . + 'Scientific' numbers are represented using + . It + uses a coefficient @c :: 'Integer'@ and a base-10 exponent @e :: 'Int'@ (do + note that since we're using an 'Int' to represent the exponent these numbers + aren't truly arbitrary precision). A scientific number corresponds to the + 'Fractional' number: @'fromInteger' c * 10 '^^' e@. + . + The main application of 'Scientific' is to be used as the target of parsing + arbitrary precision numbers coming from an untrusted source. The advantages + over using 'Rational' for this are that: + . + * 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@. + . + * 'Scientific' is safe against numbers with huge exponents. For example: + @1e1000000000 :: 'Rational'@ will fill up all space and crash your + program. Scientific works as expected: + . + > > read "1e1000000000" :: Scientific + > 1.0e1000000000 + . + * Also, the space usage of converting scientific numbers with huge exponents to + @'Integral's@ (like: 'Int') or @'RealFloat's@ (like: 'Double' or 'Float') + will always be bounded by the target type. + homepage: https://github.com/basvandijk/scientific bug-reports: https://github.com/basvandijk/scientific/issues license: BSD3 @@ -29,12 +48,17 @@ library exposed-modules: Data.Scientific + Data.Text.Lazy.Builder.Scientific + Data.ByteString.Builder.Scientific other-extensions: DeriveDataTypeable, BangPatterns ghc-options: -Wall - build-depends: base >= 4.3 && < 4.8 - , deepseq >= 1.3 && < 1.4 - , text >= 0.8 && < 1.2 - , hashable >= 1.1.2 && < 1.3 + 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 hs-source-dirs: src default-language: Haskell2010 @@ -49,8 +73,11 @@ , base >= 4.3 && < 4.8 , tasty >= 0.3.1 && < 0.9 , tasty-smallcheck >= 0.2 && < 0.9 + , tasty-quickcheck >= 0.8 && < 0.9 , smallcheck >= 1.0 && < 1.2 + , QuickCheck >= 2.7 && < 2.8 , text >= 0.8 && < 1.3 + , bytestring >= 0.10 && < 0.11 benchmark bench-scientific type: exitcode-stdio-1.0 diff --git a/src/Data/ByteString/Builder/Scientific.hs b/src/Data/ByteString/Builder/Scientific.hs new file mode 100644 --- /dev/null +++ b/src/Data/ByteString/Builder/Scientific.hs @@ -0,0 +1,134 @@ +{-# LANGUAGE CPP, MagicHash, OverloadedStrings #-} + +module Data.ByteString.Builder.Scientific + ( scientificBuilder + , formatScientificBuilder + , FPFormat(..) + ) where + +import Data.Scientific (Scientific) +import qualified Data.Scientific as Scientific + +import Data.Text.Lazy.Builder.RealFloat (FPFormat(..)) + +import qualified Data.ByteString.Char8 as BC8 + +#if !MIN_VERSION_bytestring(0,10,2) +import Data.ByteString.Lazy.Builder (Builder, string8, char8) +import Data.ByteString.Lazy.Builder.ASCII (intDec) +import Data.ByteString.Lazy.Builder.Extra (byteStringCopy) +#else +import Data.ByteString.Builder (Builder, string8, char8, intDec) +import Data.ByteString.Builder.Extra (byteStringCopy) +#endif + +import GHC.Base (Int(I#), Char(C#), chr#, ord#, (+#)) +#if MIN_VERSION_base(4,5,0) +import Data.Monoid ((<>)) +#else +import Data.Monoid (Monoid, mappend) +(<>) :: Monoid a => a -> a -> a +(<>) = mappend +infixr 6 <> +#endif + +-- | A @ByteString@ @Builder@ which renders a scientific number to full +-- precision, using standard decimal notation for arguments whose +-- absolute value lies between @0.1@ and @9,999,999@, and scientific +-- notation otherwise. +scientificBuilder :: Scientific -> Builder +scientificBuilder = formatScientificBuilder Generic Nothing + +-- | Like 'scientificBuilder' but provides rendering options. +formatScientificBuilder :: FPFormat + -> Maybe Int -- ^ Number of decimal places to render. + -> Scientific + -> Builder +formatScientificBuilder fmt decs scntfc + | scntfc < 0 = char8 '-' <> doFmt fmt (Scientific.toDecimalDigits (-scntfc)) + | otherwise = doFmt fmt (Scientific.toDecimalDigits scntfc) + where + doFmt format (is, e) = + let ds = map i2d is in + case format of + Generic -> + doFmt (if e < 0 || e > 7 then Exponent else Fixed) + (is,e) + Exponent -> + case decs of + Nothing -> + let show_e' = intDec (e-1) in + case ds of + "0" -> byteStringCopy "0.0e0" + [d] -> char8 d <> byteStringCopy ".0e" <> show_e' + (d:ds') -> char8 d <> char8 '.' <> string8 ds' <> char8 'e' <> show_e' + [] -> error $ "Data.ByteString.Builder.Scientific.formatScientificBuilder" ++ + "/doFmt/Exponent: []" + Just dec -> + let dec' = max dec 1 in + case is of + [0] -> byteStringCopy "0." <> + byteStringCopy (BC8.replicate dec' '0') <> + byteStringCopy "e0" + _ -> + let + (ei,is') = roundTo (dec'+1) is + (d:ds') = map i2d (if ei > 0 then init is' else is') + in + char8 d <> char8 '.' <> string8 ds' <> char8 'e' <> intDec (e-1+ei) + Fixed -> + let + mk0 ls = case ls of { "" -> char8 '0' ; _ -> string8 ls} + in + case decs of + Nothing + | e <= 0 -> byteStringCopy "0." <> + byteStringCopy (BC8.replicate (-e) '0') <> + string8 ds + | otherwise -> + let + f 0 s rs = mk0 (reverse s) <> char8 '.' <> mk0 rs + f n s "" = f (n-1) ('0':s) "" + f n s (r:rs) = f (n-1) (r:s) rs + in + f e "" ds + Just dec -> + let dec' = max dec 0 in + if e >= 0 then + let + (ei,is') = roundTo (dec' + e) is + (ls,rs) = splitAt (e+ei) (map i2d is') + in + mk0 ls <> (if null rs then "" else char8 '.' <> string8 rs) + else + let + (ei,is') = roundTo dec' (replicate (-e) 0 ++ is) + d:ds' = map i2d (if ei > 0 then is' else 0:is') + in + char8 d <> (if null ds' then "" else char8 '.' <> string8 ds') + +-- | Unsafe conversion for decimal digits. +{-# INLINE i2d #-} +i2d :: Int -> Char +i2d (I# i#) = C# (chr# (ord# '0'# +# i#)) + +roundTo :: Int -> [Int] -> (Int,[Int]) +roundTo d is = + case f d True is of + x@(0,_) -> x + (1,xs) -> (1, 1:xs) + _ -> error "roundTo: bad Value" + where + base = 10 + + b2 = base `quot` 2 + + f n _ [] = (0, replicate n 0) + f 0 e (x:xs) | x == b2 && e && all (== 0) xs = (0, []) -- Round to even when at exactly half the base + | otherwise = (if x >= b2 then 1 else 0, []) + f n _ (i:xs) + | i' == base = (1,0:ds) + | otherwise = (0,i':ds) + where + (c,ds) = f (n-1) (even i) xs + i' = c + i diff --git a/src/Data/Scientific.hs b/src/Data/Scientific.hs --- a/src/Data/Scientific.hs +++ b/src/Data/Scientific.hs @@ -1,7 +1,4 @@ -{-# LANGUAGE DeriveDataTypeable, BangPatterns #-} - --- TODO: The following extensions are needed for scientificBuilder: -{-# LANGUAGE CPP, MagicHash, OverloadedStrings #-} +{-# LANGUAGE DeriveDataTypeable, BangPatterns, ScopedTypeVariables #-} -- | -- Module : Data.Scientific @@ -9,62 +6,87 @@ -- License : BSD3 -- Maintainer : Bas van Dijk -- +-- @Data.Scientific@ provides a space efficient and arbitrary precision +-- scientific number type. +-- +-- 'Scientific' numbers are represented using +-- . It +-- uses an 'Integer' 'coefficient' @c@ and an 'Int' 'base10Exponent' @e@ (do +-- note that since we're using an 'Int' to represent the exponent these numbers +-- aren't truly arbitrary precision). A scientific number corresponds to the +-- 'Fractional' number: @'fromInteger' c * 10 '^^' e@. +-- +-- The main application of 'Scientific' is to be used as the target of parsing +-- arbitrary precision numbers coming from an untrusted source. The advantages +-- over using 'Rational' for this are that: +-- +-- * 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@. +-- +-- * 'Scientific' is safe against numbers with huge exponents. For example: +-- @1e1000000000 :: 'Rational'@ will fill up all space and crash your +-- program. Scientific works as expected: +-- +-- > > read "1e1000000000" :: Scientific +-- > 1.0e1000000000 +-- +-- * Also, the space usage of converting scientific numbers with huge exponents +-- to @'Integral's@ (like: 'Int') or @'RealFloat's@ (like: 'Double' or 'Float') +-- will always be bounded by the target type. +-- -- This module is designed to be imported qualified: -- -- @import Data.Scientific as Scientific@ module Data.Scientific ( Scientific + -- * Construction , scientific + -- * Projections , coefficient , base10Exponent -- * Conversions , fromFloatDigits + , toRealFloat -- * Pretty printing - , FPFormat(..) - - , scientificBuilder - , formatScientificBuilder , formatScientific + , FPFormat(..) , toDecimalDigits ) where + ---------------------------------------------------------------------- +-- Imports +---------------------------------------------------------------------- -import Control.Monad (mplus) -import Control.DeepSeq (NFData) -import Data.Char (intToDigit, ord) -import Data.Data (Data) -import Data.Function (on) -import Data.Functor ((<$>)) -import Data.Hashable (Hashable(..)) -import Data.Ratio ((%), numerator, denominator) -import Data.Typeable (Typeable) -import Numeric (floatToDigits) -import Text.Read (readPrec) +import Control.Monad (mplus) +import Control.DeepSeq (NFData) +import Data.Array (Array, listArray, (!)) +import Data.Char (intToDigit, ord) +import Data.Data (Data) +import Data.Function (on) +import Data.Functor ((<$>)) +import Data.Hashable (Hashable(..)) +import Data.Ratio ((%), numerator, denominator) +import Data.Typeable (Typeable) +import Math.NumberTheory.Logarithms (integerLog10') +import qualified Numeric (floatToDigits) +import Text.Read (readPrec) import qualified Text.ParserCombinators.ReadPrec as ReadPrec import qualified Text.ParserCombinators.ReadP as ReadP import Text.ParserCombinators.ReadP ( ReadP ) +import Data.Text.Lazy.Builder.RealFloat (FPFormat(..)) --- TODO: The following imports are needed for the scientificBuilder: -import Data.Text.Lazy.Builder (Builder, fromString, singleton, fromText) -import Data.Text.Lazy.Builder.Int (decimal) -import qualified Data.Text as T (replicate) -import GHC.Base (Int(I#), Char(C#), chr#, ord#, (+#)) -#if MIN_VERSION_base(4,5,0) -import Data.Monoid ((<>)) -#else -import Data.Monoid (Monoid, mappend) -(<>) :: Monoid a => a -> a -> a -(<>) = mappend -infixr 6 <> -#endif ---------------------------------------------------------------------- +-- Type +---------------------------------------------------------------------- -- | An arbitrary-precision number represented using -- . @@ -75,82 +97,40 @@ -- A scientific number with 'coefficient' @c@ and 'base10Exponent' @e@ -- corresponds to the 'Fractional' number: @'fromInteger' c * 10 '^^' e@ data Scientific = Scientific - { coefficient :: !Integer -- ^ The coefficient of a scientific number. - , base10Exponent :: {-# UNPACK #-} !Int -- ^ The base-10 exponent of a scientific number. + { coefficient :: !Integer + -- ^ The coefficient of a scientific number. + + , base10Exponent :: {-# UNPACK #-} !Int + -- ^ The base-10 exponent of a scientific number. } deriving (Typeable, Data) --- | @scientific c e@ constructs a scientific number with --- 'coefficient' @c@ and 'base10Exponent' @e@. +-- | @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 = 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 + + ---------------------------------------------------------------------- +-- Instances +---------------------------------------------------------------------- instance NFData Scientific instance Hashable Scientific where hashWithSalt salt = hashWithSalt salt . toRational -instance Show Scientific where - show = formatScientific Generic Nothing - -instance Read Scientific where - readPrec = ReadPrec.lift scientificP - -scientificP :: ReadP Scientific -scientificP = do - let positive = (('+' ==) <$> ReadP.satisfy isSign) `mplus` return True - pos <- positive - - let step :: Num a => a -> Int -> a - step a digit = a * 10 + fromIntegral digit - - n <- foldDigits step 0 - - let s = Scientific n 0 - fractional = foldDigits (\(Scientific a e) digit -> scientific (step a digit) (e-1)) s - - Scientific coeff expnt <- (ReadP.satisfy (== '.') >> fractional) `mplus` return s - - let signedCoeff | pos = coeff - | otherwise = negate coeff - - eP = do posE <- positive - e <- foldDigits step 0 - if posE - then return e - else return $ negate e - - (ReadP.satisfy isE >> - ((scientific signedCoeff . (expnt +)) <$> eP)) `mplus` - return (scientific signedCoeff expnt) - -foldDigits :: (a -> Int -> a) -> a -> ReadP a -foldDigits f z = ReadP.look >>= go z - where - go !a [] = return a - go !a (c:cs) - | isDecimal c = do - _ <- ReadP.get - let digit = ord c - 48 - go (f a digit) cs - | otherwise = return a - -isDecimal :: Char -> Bool -isDecimal c = c >= '0' && c <= '9' -{-# INLINE isDecimal #-} - -isSign :: Char -> Bool -isSign c = c == '-' || c == '+' -{-# INLINE isSign #-} - -isE :: Char -> Bool -isE c = c == 'e' || c == 'E' -{-# INLINE isE #-} - ----------------------------------------------------------------------- - instance Eq Scientific where (==) = (==) `on` toRational {-# INLINE (==) #-} @@ -179,40 +159,55 @@ | e1 < e2 = scientific (c1 + c2*l) e1 | otherwise = scientific (c1*r + c2 ) e2 where - l = 10 ^ (e2 - e1) - r = 10 ^ (e1 - e2) + l = magnitude (e2 - e1) + r = magnitude (e1 - e2) {-# INLINE (+) #-} Scientific c1 e1 - Scientific c2 e2 | e1 < e2 = scientific (c1 - c2*l) e1 | otherwise = scientific (c1*r - c2 ) e2 where - l = 10 ^ (e2 - e1) - r = 10 ^ (e1 - e2) + l = magnitude (e2 - e1) + r = magnitude (e1 - e2) {-# INLINE (-) #-} Scientific c1 e1 * Scientific c2 e2 = scientific (c1 * c2) (e1 + e2) {-# INLINE (*) #-} - abs (Scientific c e) = scientific (abs c) e + abs (Scientific c e) = Scientific (abs c) e {-# INLINE abs #-} - negate (Scientific c e) = scientific (negate c) e + negate (Scientific c e) = Scientific (negate c) e {-# INLINE negate #-} - signum (Scientific c _) = scientific (signum c) 0 + signum (Scientific c _) = Scientific (signum c) 0 {-# INLINE signum #-} fromInteger i = scientific i 0 {-# INLINE fromInteger #-} +-- | /WARNING:/ 'toRational' needs to compute the 'Integer' magnitude: +-- @10^e@. If applied to a huge exponent this could fill up all space +-- and crash your program! +-- +-- Avoid applying 'toRational' (or 'realToFrac') to scientific numbers +-- coming from an untrusted source and use 'toRealFloat' instead. The +-- latter guards against excessive space usage. instance Real Scientific where toRational (Scientific c e) - | e < 0 = c % (10 ^ negate e) - | otherwise = (c * 10 ^ e) % 1 + | e < 0 = c % magnitude (-e) + | otherwise = (c * magnitude e) % 1 {-# INLINE toRational #-} +{-# RULES + "realToFrac_toRealFloat_Double" + realToFrac = toRealFloat :: Scientific -> Double #-} + +{-# RULES + "realToFrac_toRealFloat_Float" + realToFrac = toRealFloat :: Scientific -> Float #-} + -- | /WARNING:/ 'recip' and '/' will diverge when their outputs have -- an infinite decimal expansion. 'fromRational' will diverge when the -- input 'Rational' has an infinite decimal expansion. @@ -220,201 +215,291 @@ recip = fromRational . recip . toRational {-# INLINE recip #-} - fromRational rational - | numer < 0 = negate $ longDiv (negate numer) 0 0 - | otherwise = longDiv numer 0 0 + fromRational rational = positivize (longDiv 0 0) (numerator rational) where - numer = numerator rational - denom = denominator rational + -- Divide the numerator by the denominator using long division. + longDiv :: Integer -> Int -> (Integer -> Scientific) + longDiv !c !e 0 = scientific c e + longDiv !c !e !n + -- TODO: Use a logarithm here! + | n < d = longDiv (c * 10) (e - 1) (n * 10) + | otherwise = longDiv (c + q) e r + where + (q, r) = n `quotRem` d - longDiv :: Integer -> Integer -> Int -> Scientific - longDiv 0 !c !e = scientific c e - longDiv !n !c !e - | n < denom = longDiv (n*10) (c * 10) (e-1) -- TODO: Use a logarithm here! - | otherwise = longDiv r (c + q) e - where - (q, r) = n `quotRem` denom + d = denominator rational instance RealFrac Scientific where - properFraction (Scientific c e) - | e < 0 = let (q, r) = c `quotRem` (10 ^ negate e) - in (fromInteger q, scientific r e) - | otherwise = (fromInteger c * 10 ^ e, 0) + -- | The function 'properFraction' takes a Scientific number @s@ + -- and returns a pair @(n,f)@ such that @s = n+f@, and: + -- + -- * @n@ is an integral number with the same sign as @s@; and + -- + -- * @f@ is a fraction with the same type and sign as @s@, + -- and with absolute value less than @1@. + properFraction s@(Scientific c e) + | e < 0 = if dangerouslySmall c e + then (0, s) + else let (q, r) = c `quotRem` magnitude (-e) + in (fromInteger q, scientific r e) + | otherwise = (fromInteger c * magnitude e, 0) {-# INLINE properFraction #-} + -- | @'truncate' s@ returns the integer nearest @s@ + -- between zero and @s@ truncate = whenFloating $ \c e -> - fromInteger $ c `quot` (10 ^ negate e) + if dangerouslySmall c e + then 0 + else fromInteger $ c `quot` magnitude (-e) {-# INLINE truncate #-} + -- | @'round' s@ returns the nearest integer to @s@; + -- the even integer if @s@ is equidistant between two integers round = whenFloating $ \c e -> - let m = c `quot` (10 ^ (negate e - 1)) - (n, r) = m `quotRem` 10 - in fromInteger $ - if c < 0 - then if r < (-5) || (r == (-5) && odd n) then n-1 else n - else if r < 5 || (r == 5 && even n) then n else n+1 + if dangerouslySmall c e + then 0 + else let (q, r) = c `quotRem` magnitude (-e) + n = fromInteger q + m = if r < 0 then n - 1 else n + 1 + f = scientific r e + in case signum $ coefficient $ abs f - 0.5 of + -1 -> n + 0 -> if even n then n else m + 1 -> m + _ -> error "round default defn: Bad value" {-# INLINE round #-} + -- | @'ceiling' s@ returns the least integer not less than @s@ ceiling = whenFloating $ \c e -> - let (q, r) = c `quotRem` (10 ^ negate e) - in fromInteger $! if r > 0 then q + 1 else q + if dangerouslySmall c e + then if c <= 0 + then 0 + else 1 + else let (q, r) = c `quotRem` magnitude (-e) + in fromInteger $! if r <= 0 then q else q + 1 {-# INLINE ceiling #-} + -- | @'floor' s@ returns the greatest integer not greater than @s@ floor = whenFloating $ \c e -> - fromInteger (c `div` (10 ^ negate e)) + if dangerouslySmall c e + then if c < 0 + then -1 + else 0 + else fromInteger (c `div` magnitude (-e)) {-# INLINE floor #-} + ---------------------------------------------------------------------- +-- Internal utilities +---------------------------------------------------------------------- +-- | This function is used in the 'RealFrac' methods to guard against +-- computing a huge magnitude (-e) which could take up all space. +-- +-- Think about parsing a scientific number from an untrusted +-- string. An attacker could supply 1e-1000000000. Lets say we want to +-- 'floor' that number to an 'Int'. When we naively try to floor it +-- using: +-- +-- @ +-- floor = whenFloating $ \c e -> +-- fromInteger (c `div` magnitude (-e)) +-- @ +-- +-- We will compute the huge Integer: @magnitude 1000000000@. This +-- computation will quickly fill up all space and crash the program. +-- +-- Note that for large /positive/ exponents there is no risk of a +-- space-leak since 'whenFloating' will compute: +-- +-- @fromInteger c * magnitude e :: a@ +-- +-- where @a@ is the target type (Int in this example). So here the +-- 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 +-- 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. +-- +-- More precisely if @dangerouslySmall c e@ returns 'True' the +-- scientific number @s@ is guaranteed to be between: +-- @-0.1 > s < 0.1@. +-- +-- Note that we avoid computing the number of decimal digits in c +-- (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 #-} + +positivize :: (Ord a, Num a, Num b) => (a -> b) -> (a -> b) +positivize f x | x < 0 = -(f (-x)) + | otherwise = f x +{-# INLINE positivize #-} + whenFloating :: (Num a) => (Integer -> Int -> a) -> Scientific -> a whenFloating f (Scientific c e) | e < 0 = f c e - | otherwise = fromInteger c * 10 ^ e + | otherwise = fromInteger c * magnitude e {-# INLINE whenFloating #-} + ---------------------------------------------------------------------- +-- Exponentiation with a cache for the most common numbers. +---------------------------------------------------------------------- --- | Efficient and exact conversion from a 'RealFloat' into a --- 'Scientific' number. -fromFloatDigits :: (RealFloat a) => a -> Scientific -fromFloatDigits rf - -- integers are way more efficient to convert via Rational. - -- We do pay the cost of always converting to Rational first though. - | denominator rat == 1 = fromInteger $ numerator rat - | rf < 0 = negate $ fromNonNegRealFloat $ negate rf - | otherwise = fromNonNegRealFloat rf +maxExpt :: Int +maxExpt = 1100 + +expts10 :: Array Int Integer +expts10 = listArray (0, maxExpt) $ iterate (*10) 1 + +-- | @magnitude e == 10 ^ e@ +magnitude :: (Num a) => Int -> a +magnitude e | e <= maxExpt = cachedPow10 e + | otherwise = cachedPow10 maxExpt * 10 ^ (e - maxExpt) where - rat = toRational rf + cachedPow10 p = fromInteger (expts10 ! p) +{-# INLINE magnitude #-} + +---------------------------------------------------------------------- +-- Conversions +---------------------------------------------------------------------- + +-- | Convert a 'RealFloat' (like a 'Double' or 'Float') into a 'Scientific' +-- number. +-- +-- Note that this function uses 'Numeric.floatToDigits' to compute the digits +-- and exponent of the 'RealFloat' number. Be aware that the algorithm used in +-- 'Numeric.floatToDigits' doesn't work as expected for some numbers, e.g. as +-- the 'Double' @1e23@ is converted to @9.9999999999999991611392e22@, and that +-- value is shown as @9.999999999999999e22@ rather than the shorter @1e23@; the +-- algorithm doesn't take the rounding direction for values exactly half-way +-- between two adjacent representable values into account, so if you have a +-- value with a short decimal representation exactly half-way between two +-- adjacent representable values, like @5^23*2^e@ for @e@ close to 23, the +-- algorithm doesn't know in which direction the short decimal representation +-- would be rounded and computes more digits +fromFloatDigits :: (RealFloat a) => a -> Scientific +fromFloatDigits = positivize fromNonNegRealFloat + where fromNonNegRealFloat r = go digits 0 0 where - (digits, e) = floatToDigits 10 r + (digits, e) = Numeric.floatToDigits 10 r - go [] !c !n = scientific c (e - n) + go [] !c !n = Scientific c (e - n) go (d:ds) !c !n = go ds (c * 10 + fromIntegral d) (n + 1) ----------------------------------------------------------------------- - --- | Similar to '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 --- --- > toDecimalDigits x = ([d1,d2,...,dn], e) --- --- then --- --- (1) @n >= 1@ +-- | Convert a 'Scientific' number into a 'RealFloat' (like a 'Double' +-- or a 'Float'). -- --- (2) @x = 0.d1d2...dn * (10^^e)@ +-- Note that this function uses 'realToFrac' +-- (@'fromRational' . 'toRational'@) internally but it guards against +-- computing huge Integer magnitudes (@10^e@) that could fill up all +-- space and crash your program. -- --- (3) @0 <= di <= 9@ -toDecimalDigits :: Scientific -> ([Int], Int) -toDecimalDigits (Scientific 0 _) = ([0], 0) -toDecimalDigits (Scientific c e) = (is, n + e) +-- Always prefer 'toRealFloat' over 'realToFrac' when converting from +-- 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 where - (is, n) = reverseAndLength $ digits c + hiLimit = ceiling (fromIntegral hi * log10Radix) + loLimit = floor (fromIntegral lo * log10Radix) - + ceiling (fromIntegral digits * log10Radix) - digits :: Integer -> [Int] - digits 0 = [] - digits i = fromIntegral r : digits q - where - (q, r) = i `quotRem` 10 + log10Radix :: Double + log10Radix = logBase 10 $ fromInteger radix - reverseAndLength :: [a] -> ([a], Int) - reverseAndLength l = rev l [] 0 - where - rev [] a !m = (a, m) - rev (x:xs) a !m = rev xs (x:a) (m+1) + radix = floatRadix (undefined :: a) + digits = floatDigits (undefined :: a) + (lo, hi) = floatRange (undefined :: a) + d = integerLog10' (abs c) + + sign x | c < 0 = -x + | otherwise = x + + ---------------------------------------------------------------------- +-- Parsing +---------------------------------------------------------------------- --- | Control the rendering of floating point numbers. -data FPFormat = Exponent - -- ^ Scientific notation (e.g. @2.3e123@). - | Fixed - -- ^ Standard decimal notation. - | Generic - -- ^ Use decimal notation for values between @0.1@ and - -- @9,999,999@, and scientific notation otherwise. - deriving (Enum, Read, Show) +instance Read Scientific where + readPrec = ReadPrec.lift scientificP --- | A @Text@ @Builder@ which renders a scientific number to full --- precision, using standard decimal notation for arguments whose --- absolute value lies between @0.1@ and @9,999,999@, and scientific --- notation otherwise. -scientificBuilder :: Scientific -> Builder -scientificBuilder = formatScientificBuilder Generic Nothing +-- A strict pair +data SP = SP !Integer {-# UNPACK #-}!Int --- | Like 'scientificBuilder' but provides rendering options. -formatScientificBuilder :: FPFormat - -> Maybe Int -- ^ Number of decimal places to render. - -> Scientific - -> Builder -formatScientificBuilder fmt decs scntfc@(Scientific c _) - | c < 0 = singleton '-' <> doFmt fmt (toDecimalDigits (-scntfc)) - | otherwise = doFmt fmt (toDecimalDigits scntfc) - where - doFmt format (is, e) = - let ds = map i2d is in - case format of - Generic -> - doFmt (if e < 0 || e > 7 then Exponent else Fixed) - (is,e) - Exponent -> - case decs of - Nothing -> - let show_e' = decimal (e-1) in - case ds of - "0" -> "0.0e0" - [d] -> singleton d <> ".0e" <> show_e' - (d:ds') -> singleton d <> singleton '.' <> fromString ds' <> singleton 'e' <> show_e' - [] -> error "formatRealFloat/doFmt/Exponent: []" - Just dec -> - let dec' = max dec 1 in - case is of - [0] -> "0." <> fromText (T.replicate dec' "0") <> "e0" - _ -> - let - (ei,is') = roundTo (dec'+1) is - (d:ds') = map i2d (if ei > 0 then init is' else is') - in - singleton d <> singleton '.' <> fromString ds' <> singleton 'e' <> decimal (e-1+ei) - Fixed -> - let - mk0 ls = case ls of { "" -> "0" ; _ -> fromString ls} - in - case decs of - Nothing - | e <= 0 -> "0." <> fromText (T.replicate (-e) "0") <> fromString ds - | otherwise -> - let - f 0 s rs = mk0 (reverse s) <> singleton '.' <> mk0 rs - f n s "" = f (n-1) ('0':s) "" - f n s (r:rs) = f (n-1) (r:s) rs - in - f e "" ds - Just dec -> - let dec' = max dec 0 in - if e >= 0 then - let - (ei,is') = roundTo (dec' + e) is - (ls,rs) = splitAt (e+ei) (map i2d is') - in - mk0 ls <> (if null rs then "" else singleton '.' <> fromString rs) - else - let - (ei,is') = roundTo dec' (replicate (-e) 0 ++ is) - d:ds' = map i2d (if ei > 0 then is' else 0:is') - in - singleton d <> (if null ds' then "" else singleton '.' <> fromString ds') +scientificP :: ReadP Scientific +scientificP = do + let positive = (('+' ==) <$> ReadP.satisfy isSign) `mplus` return True + pos <- positive --- | Unsafe conversion for decimal digits. -{-# INLINE i2d #-} -i2d :: Int -> Char -i2d (I# i#) = C# (chr# (ord# '0'# +# i#)) + let step :: Num a => a -> Int -> a + step a digit = a * 10 + fromIntegral digit + {-# INLINE step #-} + n <- foldDigits step 0 + + let s = SP n 0 + fractional = foldDigits (\(SP a e) digit -> + SP (step a digit) (e-1)) s + + SP coeff expnt <- (ReadP.satisfy (== '.') >> fractional) + `mplus` return s + + let signedCoeff | pos = coeff + | otherwise = (-coeff) + + eP = do posE <- positive + e <- foldDigits step 0 + if posE + then return e + else return (-e) + + (ReadP.satisfy isE >> + ((scientific signedCoeff . (expnt +)) <$> eP)) `mplus` + return (scientific signedCoeff expnt) + +foldDigits :: (a -> Int -> a) -> a -> ReadP a +foldDigits f z = ReadP.look >>= go z + where + go !a [] = return a + go !a (c:cs) + | isDecimal c = do + _ <- ReadP.get + let digit = ord c - 48 + go (f a digit) cs + | otherwise = return a + +isDecimal :: Char -> Bool +isDecimal c = c >= '0' && c <= '9' +{-# INLINE isDecimal #-} + +isSign :: Char -> Bool +isSign c = c == '-' || c == '+' +{-# INLINE isSign #-} + +isE :: Char -> Bool +isE c = c == 'e' || c == 'E' +{-# INLINE isE #-} + + ---------------------------------------------------------------------- +-- Pretty Printing +---------------------------------------------------------------------- +instance Show Scientific where + show = formatScientific Generic Nothing + -- | Like 'show' but provides rendering options. formatScientific :: FPFormat -> Maybe Int -- ^ Number of decimal places to render. @@ -501,3 +586,36 @@ where (c,ds) = f (n-1) (even i) xs i' = c + i + +---------------------------------------------------------------------- + +-- | Similar to '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 +-- +-- > toDecimalDigits x = ([d1,d2,...,dn], e) +-- +-- then +-- +-- (1) @n >= 1@ +-- +-- (2) @x = 0.d1d2...dn * (10^^e)@ +-- +-- (3) @0 <= di <= 9@ +toDecimalDigits :: Scientific -> ([Int], Int) +toDecimalDigits (Scientific 0 _) = ([0], 0) +toDecimalDigits (Scientific c e) = (is, n + e) + where + (is, n) = reverseAndLength $ digits c + + digits :: Integer -> [Int] + digits 0 = [] + digits i = fromIntegral r : digits q + where + (q, r) = i `quotRem` 10 + + reverseAndLength :: [a] -> ([a], Int) + reverseAndLength l = rev l [] 0 + where + rev [] a !m = (a, m) + rev (x:xs) a !m = rev xs (x:a) (m+1) diff --git a/src/Data/Text/Lazy/Builder/Scientific.hs b/src/Data/Text/Lazy/Builder/Scientific.hs new file mode 100644 --- /dev/null +++ b/src/Data/Text/Lazy/Builder/Scientific.hs @@ -0,0 +1,122 @@ +{-# LANGUAGE CPP, MagicHash, OverloadedStrings #-} + +module Data.Text.Lazy.Builder.Scientific + ( scientificBuilder + , formatScientificBuilder + , FPFormat(..) + ) where + +import Data.Scientific (Scientific) +import qualified Data.Scientific as Scientific + +import Data.Text.Lazy.Builder.RealFloat (FPFormat(..)) + +import Data.Text.Lazy.Builder (Builder, fromString, singleton, fromText) +import Data.Text.Lazy.Builder.Int (decimal) +import qualified Data.Text as T (replicate) +import GHC.Base (Int(I#), Char(C#), chr#, ord#, (+#)) +#if MIN_VERSION_base(4,5,0) +import Data.Monoid ((<>)) +#else +import Data.Monoid (Monoid, mappend) +(<>) :: Monoid a => a -> a -> a +(<>) = mappend +infixr 6 <> +#endif + +-- | A @Text@ @Builder@ which renders a scientific number to full +-- precision, using standard decimal notation for arguments whose +-- absolute value lies between @0.1@ and @9,999,999@, and scientific +-- notation otherwise. +scientificBuilder :: Scientific -> Builder +scientificBuilder = formatScientificBuilder Generic Nothing + +-- | Like 'scientificBuilder' but provides rendering options. +formatScientificBuilder :: FPFormat + -> Maybe Int -- ^ Number of decimal places to render. + -> Scientific + -> Builder +formatScientificBuilder fmt decs scntfc + | scntfc < 0 = singleton '-' <> doFmt fmt (Scientific.toDecimalDigits (-scntfc)) + | otherwise = doFmt fmt (Scientific.toDecimalDigits scntfc) + where + doFmt format (is, e) = + let ds = map i2d is in + case format of + Generic -> + doFmt (if e < 0 || e > 7 then Exponent else Fixed) + (is,e) + Exponent -> + case decs of + Nothing -> + let show_e' = decimal (e-1) in + case ds of + "0" -> "0.0e0" + [d] -> singleton d <> ".0e" <> show_e' + (d:ds') -> singleton d <> singleton '.' <> fromString ds' <> singleton 'e' <> show_e' + [] -> error $ "Data.Text.Lazy.Builder.Scientific.formatScientificBuilder" ++ + "/doFmt/Exponent: []" + Just dec -> + let dec' = max dec 1 in + case is of + [0] -> "0." <> fromText (T.replicate dec' "0") <> "e0" + _ -> + let + (ei,is') = roundTo (dec'+1) is + (d:ds') = map i2d (if ei > 0 then init is' else is') + in + singleton d <> singleton '.' <> fromString ds' <> singleton 'e' <> decimal (e-1+ei) + Fixed -> + let + mk0 ls = case ls of { "" -> "0" ; _ -> fromString ls} + in + case decs of + Nothing + | e <= 0 -> "0." <> fromText (T.replicate (-e) "0") <> fromString ds + | otherwise -> + let + f 0 s rs = mk0 (reverse s) <> singleton '.' <> mk0 rs + f n s "" = f (n-1) ('0':s) "" + f n s (r:rs) = f (n-1) (r:s) rs + in + f e "" ds + Just dec -> + let dec' = max dec 0 in + if e >= 0 then + let + (ei,is') = roundTo (dec' + e) is + (ls,rs) = splitAt (e+ei) (map i2d is') + in + mk0 ls <> (if null rs then "" else singleton '.' <> fromString rs) + else + let + (ei,is') = roundTo dec' (replicate (-e) 0 ++ is) + d:ds' = map i2d (if ei > 0 then is' else 0:is') + in + singleton d <> (if null ds' then "" else singleton '.' <> fromString ds') + +-- | Unsafe conversion for decimal digits. +{-# INLINE i2d #-} +i2d :: Int -> Char +i2d (I# i#) = C# (chr# (ord# '0'# +# i#)) + +roundTo :: Int -> [Int] -> (Int,[Int]) +roundTo d is = + case f d True is of + x@(0,_) -> x + (1,xs) -> (1, 1:xs) + _ -> error "roundTo: bad Value" + where + base = 10 + + b2 = base `quot` 2 + + f n _ [] = (0, replicate n 0) + f 0 e (x:xs) | x == b2 && e && all (== 0) xs = (0, []) -- Round to even when at exactly half the base + | otherwise = (if x >= b2 then 1 else 0, []) + f n _ (i:xs) + | i' == base = (1,0:ds) + | otherwise = (0,i':ds) + where + (c,ds) = f (n-1) (even i) xs + i' = c + i diff --git a/test/test.hs b/test/test.hs --- a/test/test.hs +++ b/test/test.hs @@ -1,29 +1,61 @@ -{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-} -{-# LANGUAGE RankNTypes, ScopedTypeVariables #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} {-# OPTIONS_GHC -fno-warn-orphans #-} module Main where -import Control.Monad -import Test.Tasty -import Test.Tasty.SmallCheck (testProperty) -import Test.SmallCheck -import Data.Scientific as Scientific -import Test.SmallCheck.Series -- (Serial, series, cons2) -import qualified Data.Text.Lazy as TL (unpack) -import qualified Data.Text.Lazy.Builder as TLB (toLazyText) +import Control.Applicative +import Control.Monad +import Data.Scientific as Scientific +import Test.Tasty +import qualified Test.SmallCheck as SC +import qualified Test.SmallCheck.Series as SC +import qualified Test.Tasty.SmallCheck as SC (testProperty) +import qualified Test.QuickCheck as QC +import qualified Test.Tasty.QuickCheck as QC (testProperty) +import qualified Data.Text.Lazy as TL (unpack) +import qualified Data.Text.Lazy.Builder as TLB (toLazyText) +import qualified Data.ByteString.Builder as B +import qualified Data.ByteString.Lazy.Char8 as BLC8 +import qualified Data.ByteString.Builder.Scientific as B +import qualified Data.Text.Lazy.Builder.Scientific as T + main :: IO () main = defaultMain $ testGroup "scientific" - [ testGroup "Formatting" + [ smallQuick "normalization" + (\s -> s /= 0 SC.==> abs (Scientific.coefficient s) `mod` 10 /= 0) + (\s -> s /= 0 QC.==> abs (Scientific.coefficient s) `mod` 10 /= 0) + + , testGroup "Formatting" [ testProperty "read . show == id" $ \s -> read (show s) === s - , testProperty "toDecimalDigits_laws" - toDecimalDigits_laws - , testProperty "Builder" $ \s -> - formatScientific Generic Nothing s == - TL.unpack (TLB.toLazyText $ formatScientificBuilder Generic Nothing s) + , smallQuick "toDecimalDigits_laws" + (SC.over nonNegativeScientificSeries toDecimalDigits_laws) + (QC.forAll nonNegativeScientificGen toDecimalDigits_laws) + + + , testGroup "Builder" + [ testProperty "Text" $ \s -> + formatScientific B.Generic Nothing s == + TL.unpack (TLB.toLazyText $ T.formatScientificBuilder B.Generic Nothing s) + + , testProperty "ByteString" $ \s -> + formatScientific B.Generic Nothing s == + BLC8.unpack (B.toLazyByteString $ B.formatScientificBuilder B.Generic Nothing s) + ] + + , testProperty "formatScientific_fromFloatDigits" $ \(d::Double) -> + formatScientific B.Generic Nothing (Scientific.fromFloatDigits d) == + show d + + -- , testProperty "formatScientific_realToFrac" $ \(d::Double) -> + -- formatScientific B.Generic Nothing (realToFrac d :: Scientific) == + -- show d ] , testGroup "Num" @@ -50,8 +82,9 @@ , testProperty "+ and negate" $ \x -> x + negate x === 0 , testProperty "- and negate" $ \x -> x - negate x === x + x - , testProperty "abs . negate == id" $ over nonNegativeScientifics $ \x -> - abs (negate x) === x + , smallQuick "abs . negate == id" + (SC.over nonNegativeScientificSeries $ \x -> abs (negate x) === x) + (QC.forAll nonNegativeScientificGen $ \x -> abs (negate x) === x) ] , testGroup "Real" @@ -88,11 +121,47 @@ ] , testGroup "Conversions" - [ testProperty "fromRealFloat" $ \s -> - Scientific.fromFloatDigits (realToFrac s :: Double) === s + [ testGroup "Float" $ conversionsProperties (undefined :: Float) + , testGroup "Double" $ conversionsProperties (undefined :: Double) ] ] +conversionsProperties :: forall realFloat. + ( RealFloat realFloat + , QC.Arbitrary realFloat + , SC.Serial IO realFloat + , Show realFloat + ) + => realFloat -> [TestTree] +conversionsProperties _ = + [ + -- testProperty "fromFloatDigits_1" $ \(d :: realFloat) -> + -- Scientific.fromFloatDigits d === realToFrac d + + -- testProperty "fromFloatDigits_2" $ \(s :: Scientific) -> + -- Scientific.fromFloatDigits (realToFrac s :: realFloat) == s + + testProperty "toRealFloat" $ \(d :: realFloat) -> + (Scientific.toRealFloat . realToFrac) d == d + + , testProperty "toRealFloat . fromFloatDigits == id" $ \(d :: realFloat) -> + (Scientific.toRealFloat . Scientific.fromFloatDigits) d == d + + -- , testProperty "fromFloatDigits . toRealFloat == id" $ \(s :: Scientific) -> + -- Scientific.fromFloatDigits (Scientific.toRealFloat s :: realFloat) == s + ] + +testProperty :: (SC.Testable IO test, QC.Testable test) + => TestName -> test -> TestTree +testProperty n test = smallQuick n test test + +smallQuick :: (SC.Testable IO smallCheck, QC.Testable quickCheck) + => TestName -> smallCheck -> quickCheck -> TestTree +smallQuick n sc qc = testGroup n + [ SC.testProperty "smallcheck" sc + , QC.testProperty "quickcheck" qc + ] + -- | ('==') specialized to 'Scientific' so we don't have to put type -- signatures everywhere. (===) :: Scientific -> Scientific -> Bool @@ -105,8 +174,8 @@ unary :: (forall a. Num a => a -> a) -> Scientific -> Bool unary op a = toRational (op a) == op (toRational a) -toDecimalDigits_laws :: (Monad m) => Property m -toDecimalDigits_laws = over nonNegativeScientifics $ \x -> +toDecimalDigits_laws :: Scientific -> Bool +toDecimalDigits_laws x = let (ds, e) = Scientific.toDecimalDigits x rule1 = length ds >= 1 @@ -152,12 +221,29 @@ _ -> error "round default defn: Bad value" ---------------------------------------------------------------------- +-- SmallCheck instances +---------------------------------------------------------------------- -instance (Monad m) => Serial m Scientific where +instance (Monad m) => SC.Serial m Scientific where series = scientifics -scientifics :: (Monad m) => Series m Scientific -scientifics = cons2 scientific +scientifics :: (Monad m) => SC.Series m Scientific +scientifics = SC.cons2 scientific -nonNegativeScientifics :: (Monad m) => Series m Scientific -nonNegativeScientifics = liftM getNonNegative series +nonNegativeScientificSeries :: (Monad m) => SC.Series m Scientific +nonNegativeScientificSeries = liftM SC.getNonNegative SC.series + + +---------------------------------------------------------------------- +-- QuickCheck instances +---------------------------------------------------------------------- + +instance QC.Arbitrary Scientific where + arbitrary = scientific <$> QC.arbitrary <*> QC.arbitrary + + shrink s = zipWith scientific (QC.shrink $ Scientific.coefficient s) + (QC.shrink $ Scientific.base10Exponent s) + +nonNegativeScientificGen :: QC.Gen Scientific +nonNegativeScientificGen = scientific <$> (QC.getNonNegative <$> QC.arbitrary) + <*> QC.arbitrary