ghc/interpreter/library/Numeric.hs

   1 #ifdef HEAD
   2 module Numeric(fromRat,
   3                showSigned, showInt,
   4                readSigned, readInt,
   5                readDec, readOct, readHex,
   6                floatToDigits,
   7                showEFloat, showFFloat, showGFloat, showFloat,
   8                readFloat, lexDigits) where
   9
  10 import Char
  11 import Array
  12
  13 import PreludeBuiltin
  14 #endif
  15 #ifdef BODY
  16
  17 -- This converts a rational to a floating.  This should be used in the
  18 -- Fractional instances of Float and Double.
  19
  20 fromRat :: (RealFloat a) => Rational -> a
  21 fromRat x =
  22     if x == 0 then encodeFloat 0 0              -- Handle exceptional cases
  23     else if x < 0 then - fromRat' (-x)          -- first.
  24     else fromRat' x
  25
  26 -- Conversion process:
  27 -- Scale the rational number by the RealFloat base until
  28 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
  29 -- Then round the rational to an Integer and encode it with the exponent
  30 -- that we got from the scaling.
  31 -- To speed up the scaling process we compute the log2 of the number to get
  32 -- a first guess of the exponent.
  33 fromRat' :: (RealFloat a) => Rational -> a
  34 fromRat' x = r
  35   where b = floatRadix r
  36         p = floatDigits r
  37         (minExp0, _) = floatRange r
  38         minExp = minExp0 - p            -- the real minimum exponent
  39         xMin = toRational (expt b (p-1))
  40         xMax = toRational (expt b p)
  41         p0 = (integerLogBase b (numerator x) -
  42               integerLogBase b (denominator x) - p) `max` minExp
  43         f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
  44         (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
  45         r = encodeFloat (round x') p'
  46
  47 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
  48 scaleRat :: Rational -> Int -> Rational -> Rational ->
  49              Int -> Rational -> (Rational, Int)
  50 scaleRat b minExp xMin xMax p x
  51     | p <= minExp = (x, p)
  52     | x >= xMax   = scaleRat b minExp xMin xMax (p+1) (x/b)
  53     | x < xMin    = scaleRat b minExp xMin xMax (p-1) (x*b)
  54     | otherwise   = (x, p)
  55
  56 -- Exponentiation with a cache for the most common numbers.
  57 minExpt = 0::Int
  58 maxExpt = 1100::Int
  59 expt :: BIGNUMTYPE -> Int -> BIGNUMTYPE
  60 expt base n =
  61     if base == 2 && n >= minExpt && n <= maxExpt then
  62         expts!n
  63     else
  64         base^n
  65
  66 expts :: Array Int BIGNUMTYPE
  67 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
  68
  69 -- Compute the (floor of the) log of i in base b.
  70 -- Simplest way would be just divide i by b until it's smaller then b,
  71 -- but that would be very slow!  We are just slightly more clever.
  72 integerLogBase :: BIGNUMTYPE -> BIGNUMTYPE -> Int
  73 integerLogBase b i =
  74      if i < b then
  75         0
  76      else
  77         -- Try squaring the base first to cut down the number of divisions.
  78         let l = 2 * integerLogBase (b*b) i
  79             doDiv :: BIGNUMTYPE -> Int -> Int
  80             doDiv i l = if i < b then l else doDiv (i `div` b) (l+1)
  81         in  doDiv (i `div` (b^l)) l
  82
  83
  84 -- Misc utilities to show integers and floats
  85
  86 showSigned    :: Real a => (a -> ShowS) -> Int -> a -> ShowS
  87 showSigned showPos p x | x < 0 = showParen (p > 6)
  88                                            (showChar '-' . showPos (-x))
  89                        | otherwise = showPos x
  90
  91 -- showInt is used for positive numbers only
  92 showInt    :: Integral a => a -> ShowS
  93 showInt n r | n < 0 = error "Numeric.showInt: can't show negative numbers"
  94             | otherwise =
  95               let (n',d) = quotRem n 10
  96                   r'     = toEnum (fromEnum '0' + fromIntegral d) : r
  97               in  if n' == 0 then r' else showInt n' r'
  98
  99
 100 readSigned :: (Real a) => ReadS a -> ReadS a
 101 readSigned readPos = readParen False read'
 102                      where read' r  = read'' r ++
 103                                       [(-x,t) | ("-",s) <- lex r,
 104                                                 (x,t)   <- read'' s]
 105                            read'' r = [(n,s)  | (str,s) <- lex r,
 106                                                 (n,"")  <- readPos str]
 107
 108
 109 -- readInt reads a string of digits using an arbitrary base.
 110 -- Leading minus signs must be handled elsewhere.
 111
 112 readInt :: (Integral a) => a -> (Char -> Bool) -> (Char -> Int) -> ReadS a
 113 readInt radix isDig digToInt s =
 114    [(foldl1 (\n d -> n * radix + d) (map (fromIntegral . digToInt) ds), r)
 115           | (ds,r) <- nonnull isDig s ]
 116
 117 -- Unsigned readers for various bases
 118 readDec, readOct, readHex :: (Integral a) => ReadS a
 119 readDec = readInt 10 isDigit digitToInt
 120 readOct = readInt  8 isOctDigit digitToInt
 121 readHex = readInt 16 isHexDigit digitToInt
 122
 123
 124 showEFloat     :: (RealFloat a) => Maybe Int -> a -> ShowS
 125 showFFloat     :: (RealFloat a) => Maybe Int -> a -> ShowS
 126 showGFloat     :: (RealFloat a) => Maybe Int -> a -> ShowS
 127 showFloat      :: (RealFloat a) => a -> ShowS
 128
 129 showEFloat d x =  showString (formatRealFloat FFExponent d x)
 130 showFFloat d x =  showString (formatRealFloat FFFixed d x)
 131 showGFloat d x =  showString (formatRealFloat FFGeneric d x)
 132 showFloat      =  showGFloat Nothing
 133
 134 -- These are the format types.  This type is not exported.
 135
 136 data FFFormat = FFExponent | FFFixed | FFGeneric
 137
 138 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
 139 formatRealFloat fmt decs x = s
 140   where base = 10
 141         s = if isNaN x then
 142                 "NaN"
 143             else if isInfinite x then
 144                 if x < 0 then "-Infinity" else "Infinity"
 145             else if x < 0 || isNegativeZero x then
 146                 '-' : doFmt fmt (floatToDigits (toInteger base) (-x))
 147             else
 148                 doFmt fmt (floatToDigits (toInteger base) x)
 149         doFmt fmt (is, e) =
 150             let ds = map intToDigit is
 151             in  case fmt of
 152                 FFGeneric ->
 153                     doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
 154                           (is, e)
 155                 FFExponent ->
 156                     case decs of
 157                     Nothing ->
 158                         case ds of
 159                          ['0'] -> "0.0e0"
 160                          [d]   -> d : ".0e" ++ show (e-1)
 161                          d:ds  -> d : '.' : ds ++ 'e':show (e-1)
 162                     Just dec ->
 163                         let dec' = max dec 1 in
 164                         case is of
 165                          [0] -> '0':'.':take dec' (repeat '0') ++ "e0"
 166                          _ ->
 167                           let (ei, is') = roundTo base (dec'+1) is
 168                               d:ds = map intToDigit
 169                                          (if ei > 0 then init is' else is')
 170                           in d:'.':ds  ++ "e" ++ show (e-1+ei)
 171                 FFFixed ->
 172                     case decs of
 173                     Nothing ->
 174                         let f 0 s ds = mk0 s ++ "." ++ mk0 ds
 175                             f n s "" = f (n-1) (s++"0") ""
 176                             f n s (d:ds) = f (n-1) (s++[d]) ds
 177                             mk0 "" = "0"
 178                             mk0 s = s
 179                         in  f e "" ds
 180                     Just dec ->
 181                         let dec' = max dec 0 in
 182                         if e >= 0 then
 183                             let (ei, is') = roundTo base (dec' + e) is
 184                                 (ls, rs) = splitAt (e+ei) (map intToDigit is')
 185                             in  (if null ls then "0" else ls) ++
 186                                 (if null rs then "" else '.' : rs)
 187                         else
 188                             let (ei, is') = roundTo base dec'
 189                                               (replicate (-e) 0 ++ is)
 190                                 d : ds = map intToDigit
 191                                             (if ei > 0 then is' else 0:is')
 192                             in  d : '.' : ds
 193
 194 roundTo :: Int -> Int -> [Int] -> (Int, [Int])
 195 roundTo base d is = case f d is of
 196                 (0, is) -> (0, is)
 197                 (1, is) -> (1, 1 : is)
 198   where b2 = base `div` 2
 199         f n [] = (0, replicate n 0)
 200         f 0 (i:_) = (if i >= b2 then 1 else 0, [])
 201         f d (i:is) =
 202             let (c, ds) = f (d-1) is
 203                 i' = c + i
 204             in  if i' == base then (1, 0:ds) else (0, i':ds)
 205
 206 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
 207 -- by R.G. Burger and R. K. Dybvig, in PLDI 96.
 208 -- This version uses a much slower logarithm estimator.  It should be improved.
 209
 210 -- This function returns a list of digits (Ints in [0..base-1]) and an
 211 -- exponent.
 212
 213 floatToDigits :: (RealFloat a) => BIGNUMTYPE -> a -> ([Int], Int)
 214
 215 floatToDigits _ 0 = ([0], 0)
 216 floatToDigits base x =
 217     let (f0, e0) = decodeFloat x
 218         (minExp0, _) = floatRange x
 219         p = floatDigits x
 220         b = floatRadix x
 221         minExp = minExp0 - p            -- the real minimum exponent
 222         -- Haskell requires that f be adjusted so denormalized numbers
 223         -- will have an impossibly low exponent.  Adjust for this.
 224         (f, e) = let n = minExp - e0
 225                  in  if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
 226
 227         (r, s, mUp, mDn) =
 228            if e >= 0 then
 229                let be = b^e in
 230                if f == b^(p-1) then
 231                    (f*be*b*2, 2*b, be*b, b)
 232                else
 233                    (f*be*2, 2, be, be)
 234            else
 235                if e > minExp && f == b^(p-1) then
 236                    (f*b*2, b^(-e+1)*2, b, 1)
 237                else
 238                    (f*2, b^(-e)*2, 1, 1)
 239         k =
 240             let k0 =
 241 #if 1 /* hack to overcome temporary Hugs bug (fixed size Integers) */
 242                      0
 243 #else
 244                     if b==2 && base==10 then
 245                         -- logBase 10 2 is slightly bigger than 3/10 so
 246                         -- the following will err on the low side.  Ignoring
 247                         -- the fraction will make it err even more.
 248                         -- Haskell promises that p-1 <= logBase b f < p.
 249                         (p - 1 + e0) * 3 `div` 10
 250                     else
 251                         ceiling ((log (fromInteger (f+1)) +
 252                                  fromInt e * log (fromInteger b)) /
 253                                   log (fromInteger base) `asTypeOf` x)
 254 #endif
 255                 fixup n =
 256                     if n >= 0 then
 257                         if r + mUp <= expt base n * s then n else fixup (n+1)
 258                     else
 259                         if expt base (-n) * (r + mUp) <= s then n
 260                                                            else fixup (n+1)
 261             in  fixup k0
 262
 263         gen ds rn sN mUpN mDnN =
 264             let (dn, rn') = (rn * base) `divMod` sN
 265                 mUpN' = mUpN * base
 266                 mDnN' = mDnN * base
 267             in  case (rn' < mDnN', rn' + mUpN' > sN) of
 268                 (True,  False) -> dn : ds
 269                 (False, True)  -> dn+1 : ds
 270                 (True,  True)  -> if rn' * 2 < sN then dn : ds else dn+1 : ds
 271                 (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
 272         rds =
 273             if k >= 0 then
 274                 gen [] r (s * expt base k) mUp mDn
 275             else
 276                 let bk = expt base (-k)
 277                 in  gen [] (r * bk) s (mUp * bk) (mDn * bk)
 278     in  (map toInt (reverse rds), k)
 279
 280
 281
 282 -- This floating point reader uses a less restrictive syntax for floating
 283 -- point than the Haskell lexer.  The `.' is optional.
 284
 285 readFloat     :: (RealFloat a) => ReadS a
 286 readFloat r    = [(fromRational ((n%1)*10^^(k-d)),t) | (n,d,s) <- readFix r,
 287                                                        (k,t)   <- readExp s]
 288                  where readFix r = [(read (ds++ds'), length ds', t)
 289                                         | (ds,d) <- lexDigits r,
 290                                           (ds',t) <- lexFrac d ]
 291
 292                        lexFrac ('.':ds) = lexDigits ds
 293                        lexFrac s        = [("",s)]
 294
 295                        readExp (e:s) | e `elem` "eE" = readExp' s
 296                        readExp s                     = [(0,s)]
 297
 298                        readExp' ('-':s) = [(-k,t) | (k,t) <- readDec s]
 299                        readExp' ('+':s) = readDec s
 300                        readExp' s       = readDec s
 301
 302 lexDigits        :: ReadS String
 303 lexDigits        =  nonnull isDigit
 304
 305 nonnull          :: (Char -> Bool) -> ReadS String
 306 nonnull p s      =  [(cs,t) | (cs@(_:_),t) <- [span p s]]
 307
 308 #endif /* BODY */