X-Git-Url: http://git.megacz.com/?a=blobdiff_plain;f=Numeric.hs;h=4b202d0efd67d0705eceaaba1f45dde6eeb5caf6;hb=7dbb606d7b57cdad87a0ffbdb6ea4a274ebca7c0;hp=01036b22bc949cc6d9479f67fa19d0f24f9b7b05;hpb=9fa9bc17072a58c0bae2cce4764d38677e96ac29;p=ghc-base.git diff --git a/Numeric.hs b/Numeric.hs index 01036b2..4b202d0 100644 --- a/Numeric.hs +++ b/Numeric.hs @@ -1,51 +1,61 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE CPP, NoImplicitPrelude, MagicHash #-} + ----------------------------------------------------------------------------- -- | -- Module : Numeric -- Copyright : (c) The University of Glasgow 2002 --- License : BSD-style (see the file libraries/core/LICENSE) +-- License : BSD-style (see the file libraries/base/LICENSE) -- -- Maintainer : libraries@haskell.org -- Stability : provisional -- Portability : portable -- --- $Id: Numeric.hs,v 1.7 2002/04/24 16:31:37 simonmar Exp $ --- -- Odds and ends, mostly functions for reading and showing --- RealFloat-like kind of values. +-- 'RealFloat'-like kind of values. -- ----------------------------------------------------------------------------- module Numeric ( - fromRat, -- :: (RealFloat a) => Rational -> a - showSigned, -- :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS - readSigned, -- :: (Real a) => ReadS a -> ReadS a + -- * Showing - readInt, -- :: (Integral a) => a -> (Char -> Bool) - -- -> (Char -> Int) -> ReadS a - readDec, -- :: (Integral a) => ReadS a - readOct, -- :: (Integral a) => ReadS a - readHex, -- :: (Integral a) => ReadS a + showSigned, -- :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS - showInt, -- :: Integral a => a -> ShowS showIntAtBase, -- :: Integral a => a -> (a -> Char) -> a -> ShowS + showInt, -- :: Integral a => a -> ShowS showHex, -- :: Integral a => a -> ShowS showOct, -- :: Integral a => a -> ShowS - showBin, -- :: Integral a => a -> ShowS - showEFloat, -- :: (RealFloat a) => Maybe Int -> a -> ShowS - showFFloat, -- :: (RealFloat a) => Maybe Int -> a -> ShowS - showGFloat, -- :: (RealFloat a) => Maybe Int -> a -> ShowS - showFloat, -- :: (RealFloat a) => a -> ShowS - readFloat, -- :: (RealFloat a) => ReadS a - - floatToDigits, -- :: (RealFloat a) => Integer -> a -> ([Int], Int) - lexDigits, -- :: ReadS String + showEFloat, -- :: (RealFloat a) => Maybe Int -> a -> ShowS + showFFloat, -- :: (RealFloat a) => Maybe Int -> a -> ShowS + showGFloat, -- :: (RealFloat a) => Maybe Int -> a -> ShowS + showFloat, -- :: (RealFloat a) => a -> ShowS - ) where + floatToDigits, -- :: (RealFloat a) => Integer -> a -> ([Int], Int) -import Data.Char + -- * Reading + + -- | /NB:/ 'readInt' is the \'dual\' of 'showIntAtBase', + -- and 'readDec' is the \`dual\' of 'showInt'. + -- The inconsistent naming is a historical accident. + + readSigned, -- :: (Real a) => ReadS a -> ReadS a + + readInt, -- :: (Integral a) => a -> (Char -> Bool) + -- -> (Char -> Int) -> ReadS a + readDec, -- :: (Integral a) => ReadS a + readOct, -- :: (Integral a) => ReadS a + readHex, -- :: (Integral a) => ReadS a + + readFloat, -- :: (RealFloat a) => ReadS a + + lexDigits, -- :: ReadS String + + -- * Miscellaneous + + fromRat, -- :: (RealFloat a) => Rational -> a + + ) where #ifdef __GLASGOW_HASKELL__ import GHC.Base @@ -57,66 +67,78 @@ import GHC.Show import Data.Maybe import Text.ParserCombinators.ReadP( ReadP, readP_to_S, pfail ) import qualified Text.Read.Lex as L +#else +import Data.Char #endif #ifdef __HUGS__ -import Array +import Hugs.Prelude +import Hugs.Numeric #endif - --- ********************************************************* --- * * --- \subsection{Reading} --- * * --- ********************************************************* - -readInt :: Num a => a -> (Char -> Bool) -> (Char -> Int) -> ReadS a +#ifdef __GLASGOW_HASKELL__ +-- ----------------------------------------------------------------------------- +-- Reading + +-- | Reads an /unsigned/ 'Integral' value in an arbitrary base. +readInt :: Num a + => a -- ^ the base + -> (Char -> Bool) -- ^ a predicate distinguishing valid digits in this base + -> (Char -> Int) -- ^ a function converting a valid digit character to an 'Int' + -> ReadS a readInt base isDigit valDigit = readP_to_S (L.readIntP base isDigit valDigit) -readOct, readDec, readHex :: Num a => ReadS a +-- | Read an unsigned number in octal notation. +readOct :: Num a => ReadS a readOct = readP_to_S L.readOctP + +-- | Read an unsigned number in decimal notation. +readDec :: Num a => ReadS a readDec = readP_to_S L.readDecP + +-- | Read an unsigned number in hexadecimal notation. +-- Both upper or lower case letters are allowed. +readHex :: Num a => ReadS a readHex = readP_to_S L.readHexP +-- | Reads an /unsigned/ 'RealFrac' value, +-- expressed in decimal scientific notation. readFloat :: RealFrac a => ReadS a readFloat = readP_to_S readFloatP readFloatP :: RealFrac a => ReadP a readFloatP = - do L.Number x <- L.lex - case L.numberToRational x of - Nothing -> pfail - Just y -> return (fromRational y) + do tok <- L.lex + case tok of + L.Rat y -> return (fromRational y) + L.Int i -> return (fromInteger i) + _ -> pfail -- It's turgid to have readSigned work using list comprehensions, -- but it's specified as a ReadS to ReadS transformer -- With a bit of luck no one will use it. + +-- | Reads a /signed/ 'Real' value, given a reader for an unsigned value. readSigned :: (Real a) => ReadS a -> ReadS a readSigned readPos = readParen False read' - where read' r = read'' r ++ - (do - ("-",s) <- lex r - (x,t) <- read'' s - return (-x,t)) - read'' r = do - (str,s) <- lex r - (n,"") <- readPos str - return (n,s) - - --- ********************************************************* --- * * --- \subsection{Showing} --- * * --- ********************************************************* - - - -#ifdef __GLASGOW_HASKELL__ + where read' r = read'' r ++ + (do + ("-",s) <- lex r + (x,t) <- read'' s + return (-x,t)) + read'' r = do + (str,s) <- lex r + (n,"") <- readPos str + return (n,s) + +-- ----------------------------------------------------------------------------- +-- Showing + +-- | Show /non-negative/ 'Integral' numbers in base 10. showInt :: Integral a => a -> ShowS -showInt n cs - | n < 0 = error "Numeric.showInt: can't show negative numbers" - | otherwise = go n cs +showInt n0 cs0 + | n0 < 0 = error "Numeric.showInt: can't show negative numbers" + | otherwise = go n0 cs0 where go n cs | n < 10 = case unsafeChr (ord '0' + fromIntegral n) of @@ -131,285 +153,68 @@ showInt n cs -- mutual module deps. {-# SPECIALIZE showEFloat :: - Maybe Int -> Float -> ShowS, - Maybe Int -> Double -> ShowS #-} + Maybe Int -> Float -> ShowS, + Maybe Int -> Double -> ShowS #-} {-# SPECIALIZE showFFloat :: - Maybe Int -> Float -> ShowS, - Maybe Int -> Double -> ShowS #-} + Maybe Int -> Float -> ShowS, + Maybe Int -> Double -> ShowS #-} {-# SPECIALIZE showGFloat :: - Maybe Int -> Float -> ShowS, - Maybe Int -> Double -> ShowS #-} + Maybe Int -> Float -> ShowS, + Maybe Int -> Double -> ShowS #-} +-- | Show a signed 'RealFloat' value +-- using scientific (exponential) notation (e.g. @2.45e2@, @1.5e-3@). +-- +-- In the call @'showEFloat' digs val@, if @digs@ is 'Nothing', +-- the value is shown to full precision; if @digs@ is @'Just' d@, +-- then at most @d@ digits after the decimal point are shown. showEFloat :: (RealFloat a) => Maybe Int -> a -> ShowS -showFFloat :: (RealFloat a) => Maybe Int -> a -> ShowS -showGFloat :: (RealFloat a) => Maybe Int -> a -> ShowS -showEFloat d x = showString (formatRealFloat FFExponent d x) -showFFloat d x = showString (formatRealFloat FFFixed d x) -showGFloat d x = showString (formatRealFloat FFGeneric d x) -#endif +-- | Show a signed 'RealFloat' value +-- using standard decimal notation (e.g. @245000@, @0.0015@). +-- +-- In the call @'showFFloat' digs val@, if @digs@ is 'Nothing', +-- the value is shown to full precision; if @digs@ is @'Just' d@, +-- then at most @d@ digits after the decimal point are shown. +showFFloat :: (RealFloat a) => Maybe Int -> a -> ShowS -#ifdef __HUGS__ --- This converts a rational to a floating. This should be used in the --- Fractional instances of Float and Double. - -fromRat :: (RealFloat a) => Rational -> a -fromRat x = - if x == 0 then encodeFloat 0 0 -- Handle exceptional cases - else if x < 0 then - fromRat' (-x) -- first. - else fromRat' x - --- Conversion process: --- Scale the rational number by the RealFloat base until --- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat). --- Then round the rational to an Integer and encode it with the exponent --- that we got from the scaling. --- To speed up the scaling process we compute the log2 of the number to get --- a first guess of the exponent. -fromRat' :: (RealFloat a) => Rational -> a -fromRat' x = r - where b = floatRadix r - p = floatDigits r - (minExp0, _) = floatRange r - minExp = minExp0 - p -- the real minimum exponent - xMin = toRational (expt b (p-1)) - xMax = toRational (expt b p) - p0 = (integerLogBase b (numerator x) - - integerLogBase b (denominator x) - p) `max` minExp - f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1 - (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f) - r = encodeFloat (round x') p' - --- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp. -scaleRat :: Rational -> Int -> Rational -> Rational -> - Int -> Rational -> (Rational, Int) -scaleRat b minExp xMin xMax p x = - if p <= minExp then - (x, p) - else if x >= xMax then - scaleRat b minExp xMin xMax (p+1) (x/b) - else if x < xMin then - scaleRat b minExp xMin xMax (p-1) (x*b) - else - (x, p) - --- Exponentiation with a cache for the most common numbers. -minExpt = 0::Int -maxExpt = 1100::Int -expt :: Integer -> Int -> Integer -expt base n = - if base == 2 && n >= minExpt && n <= maxExpt then - expts!n - else - base^n - -expts :: Array Int Integer -expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]] - --- Compute the (floor of the) log of i in base b. --- Simplest way would be just divide i by b until it's smaller then b, --- but that would be very slow! We are just slightly more clever. -integerLogBase :: Integer -> Integer -> Int -integerLogBase b i = - if i < b then - 0 - else - -- Try squaring the base first to cut down the number of divisions. - let l = 2 * integerLogBase (b*b) i - doDiv :: Integer -> Int -> Int - doDiv i l = if i < b then l else doDiv (i `div` b) (l+1) - in doDiv (i `div` (b^l)) l - - --- Misc utilities to show integers and floats - -showEFloat :: (RealFloat a) => Maybe Int -> a -> ShowS -showFFloat :: (RealFloat a) => Maybe Int -> a -> ShowS -showGFloat :: (RealFloat a) => Maybe Int -> a -> ShowS -showFloat :: (RealFloat a) => a -> ShowS +-- | Show a signed 'RealFloat' value +-- using standard decimal notation for arguments whose absolute value lies +-- between @0.1@ and @9,999,999@, and scientific notation otherwise. +-- +-- In the call @'showGFloat' digs val@, if @digs@ is 'Nothing', +-- the value is shown to full precision; if @digs@ is @'Just' d@, +-- then at most @d@ digits after the decimal point are shown. +showGFloat :: (RealFloat a) => Maybe Int -> a -> ShowS showEFloat d x = showString (formatRealFloat FFExponent d x) showFFloat d x = showString (formatRealFloat FFFixed d x) showGFloat d x = showString (formatRealFloat FFGeneric d x) -showFloat = showGFloat Nothing - --- These are the format types. This type is not exported. - -data FFFormat = FFExponent | FFFixed | FFGeneric - -formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String -formatRealFloat fmt decs x = s - where base = 10 - s = if isNaN x then - "NaN" - else if isInfinite x then - if x < 0 then "-Infinity" else "Infinity" - else if x < 0 || isNegativeZero x then - '-' : doFmt fmt (floatToDigits (toInteger base) (-x)) - else - doFmt fmt (floatToDigits (toInteger base) x) - doFmt fmt (is, e) = - let ds = map intToDigit is - in case fmt of - FFGeneric -> - doFmt (if e < 0 || e > 7 then FFExponent else FFFixed) - (is, e) - FFExponent -> - case decs of - Nothing -> - case ds of - ['0'] -> "0.0e0" - [d] -> d : ".0e" ++ show (e-1) - d:ds -> d : '.' : ds ++ 'e':show (e-1) - Just dec -> - let dec' = max dec 1 in - case is of - [0] -> '0':'.':take dec' (repeat '0') ++ "e0" - _ -> - let (ei, is') = roundTo base (dec'+1) is - d:ds = map intToDigit - (if ei > 0 then init is' else is') - in d:'.':ds ++ "e" ++ show (e-1+ei) - FFFixed -> - case decs of - Nothing -> - let f 0 s ds = mk0 s ++ "." ++ mk0 ds - f n s "" = f (n-1) (s++"0") "" - f n s (d:ds) = f (n-1) (s++[d]) ds - mk0 "" = "0" - mk0 s = s - in f e "" ds - Just dec -> - let dec' = max dec 0 in - if e >= 0 then - let (ei, is') = roundTo base (dec' + e) is - (ls, rs) = splitAt (e+ei) (map intToDigit is') - in (if null ls then "0" else ls) ++ - (if null rs then "" else '.' : rs) - else - let (ei, is') = roundTo base dec' - (replicate (-e) 0 ++ is) - d : ds = map intToDigit - (if ei > 0 then is' else 0:is') - in d : '.' : ds - -roundTo :: Int -> Int -> [Int] -> (Int, [Int]) -roundTo base d is = case f d is of - (0, is) -> (0, is) - (1, is) -> (1, 1 : is) - where b2 = base `div` 2 - f n [] = (0, replicate n 0) - f 0 (i:_) = (if i >= b2 then 1 else 0, []) - f d (i:is) = - let (c, ds) = f (d-1) is - i' = c + i - in if i' == base then (1, 0:ds) else (0, i':ds) - --- --- Based on "Printing Floating-Point Numbers Quickly and Accurately" --- by R.G. Burger and R. K. Dybvig, in PLDI 96. --- This version uses a much slower logarithm estimator. It should be improved. - --- This function returns a list of digits (Ints in [0..base-1]) and an --- exponent. - -floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int) - -floatToDigits _ 0 = ([0], 0) -floatToDigits base x = - let (f0, e0) = decodeFloat x - (minExp0, _) = floatRange x - p = floatDigits x - b = floatRadix x - minExp = minExp0 - p -- the real minimum exponent - -- Haskell requires that f be adjusted so denormalized numbers - -- will have an impossibly low exponent. Adjust for this. - (f, e) = let n = minExp - e0 - in if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0) - - (r, s, mUp, mDn) = - if e >= 0 then - let be = b^e in - if f == b^(p-1) then - (f*be*b*2, 2*b, be*b, b) - else - (f*be*2, 2, be, be) - else - if e > minExp && f == b^(p-1) then - (f*b*2, b^(-e+1)*2, b, 1) - else - (f*2, b^(-e)*2, 1, 1) - k = - let k0 = - if b==2 && base==10 then - -- logBase 10 2 is slightly bigger than 3/10 so - -- the following will err on the low side. Ignoring - -- the fraction will make it err even more. - -- Haskell promises that p-1 <= logBase b f < p. - (p - 1 + e0) * 3 `div` 10 - else - ceiling ((log (fromInteger (f+1)) + - fromIntegral e * log (fromInteger b)) / - log (fromInteger base)) - fixup n = - if n >= 0 then - if r + mUp <= expt base n * s then n else fixup (n+1) - else - if expt base (-n) * (r + mUp) <= s then n - else fixup (n+1) - in fixup k0 - - gen ds rn sN mUpN mDnN = - let (dn, rn') = (rn * base) `divMod` sN - mUpN' = mUpN * base - mDnN' = mDnN * base - in case (rn' < mDnN', rn' + mUpN' > sN) of - (True, False) -> dn : ds - (False, True) -> dn+1 : ds - (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds - (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN' - rds = - if k >= 0 then - gen [] r (s * expt base k) mUp mDn - else - let bk = expt base (-k) - in gen [] (r * bk) s (mUp * bk) (mDn * bk) - in (map fromIntegral (reverse rds), k) -#endif +#endif /* __GLASGOW_HASKELL__ */ -- --------------------------------------------------------------------------- -- Integer printing functions -showIntAtBase :: Integral a => a -> (a -> Char) -> a -> ShowS -showIntAtBase base toChr n r - | n < 0 = error ("Numeric.showIntAtBase: applied to negative number " ++ show n) - | otherwise = - case quotRem n base of { (n', d) -> - let c = toChr d in - seq c $ -- stricter than necessary - let - r' = c : r - in - if n' == 0 then r' else showIntAtBase base toChr n' r' - } - +-- | Shows a /non-negative/ 'Integral' number using the base specified by the +-- first argument, and the character representation specified by the second. +showIntAtBase :: Integral a => a -> (Int -> Char) -> a -> ShowS +showIntAtBase base toChr n0 r0 + | base <= 1 = error ("Numeric.showIntAtBase: applied to unsupported base " ++ show base) + | n0 < 0 = error ("Numeric.showIntAtBase: applied to negative number " ++ show n0) + | otherwise = showIt (quotRem n0 base) r0 + where + showIt (n,d) r = seq c $ -- stricter than necessary + case n of + 0 -> r' + _ -> showIt (quotRem n base) r' + where + c = toChr (fromIntegral d) + r' = c : r + +-- | Show /non-negative/ 'Integral' numbers in base 16. showHex :: Integral a => a -> ShowS -showHex n r = - showString "0x" $ - showIntAtBase 16 (toChrHex) n r - where - toChrHex d - | d < 10 = chr (ord '0' + fromIntegral d) - | otherwise = chr (ord 'a' + fromIntegral (d - 10)) +showHex = showIntAtBase 16 intToDigit +-- | Show /non-negative/ 'Integral' numbers in base 8. showOct :: Integral a => a -> ShowS -showOct n r = - showString "0o" $ - showIntAtBase 8 (toChrOct) n r - where toChrOct d = chr (ord '0' + fromIntegral d) - -showBin :: Integral a => a -> ShowS -showBin n r = - showString "0b" $ - showIntAtBase 2 (toChrOct) n r - where toChrOct d = chr (ord '0' + fromIntegral d) +showOct = showIntAtBase 8 intToDigit