ghc/lib/std/PrelNumExtra.lhs

   1 %
   2 % (c) The AQUA Project, Glasgow University, 1994-1996
   3 %
   4
   5 \section[PrelNumExtra]{Module @PrelNumExtra@}
   6
   7 \begin{code}
   8 {-# OPTIONS -fno-implicit-prelude #-}
   9 {-# OPTIONS -H20m #-}
  10
  11 #include "../includes/ieee-flpt.h"
  12
  13 \end{code}
  14
  15 \begin{code}
  16 module PrelNumExtra where
  17
  18 import PrelBase
  19 import PrelGHC
  20 import PrelEnum
  21 import PrelShow
  22 import PrelNum
  23 import PrelErr ( error )
  24 import PrelList
  25 import PrelMaybe
  26 import Maybe            ( fromMaybe )
  27
  28 import PrelArr          ( Array, array, (!) )
  29 import PrelIOBase       ( unsafePerformIO )
  30 import PrelCCall        ()      -- we need the definitions of CCallable and
  31                                 -- CReturnable for the _ccall_s herein.
  32 \end{code}
  33
  34 %*********************************************************
  35 %*                                                      *
  36 \subsection{Type @Float@}
  37 %*                                                      *
  38 %*********************************************************
  39
  40 \begin{code}
  41 instance Eq Float where
  42     (F# x) == (F# y) = x `eqFloat#` y
  43
  44 instance Ord Float where
  45     (F# x) `compare` (F# y) | x `ltFloat#` y = LT
  46                             | x `eqFloat#` y = EQ
  47                             | otherwise      = GT
  48
  49     (F# x) <  (F# y) = x `ltFloat#`  y
  50     (F# x) <= (F# y) = x `leFloat#`  y
  51     (F# x) >= (F# y) = x `geFloat#`  y
  52     (F# x) >  (F# y) = x `gtFloat#`  y
  53
  54 instance  Num Float  where
  55     (+)         x y     =  plusFloat x y
  56     (-)         x y     =  minusFloat x y
  57     negate      x       =  negateFloat x
  58     (*)         x y     =  timesFloat x y
  59     abs x | x >= 0.0    =  x
  60           | otherwise   =  negateFloat x
  61     signum x | x == 0.0  = 0
  62              | x > 0.0   = 1
  63              | otherwise = negate 1
  64
  65     {-# INLINE fromInteger #-}
  66     fromInteger n       =  encodeFloat n 0
  67         -- It's important that encodeFloat inlines here, and that
  68         -- fromInteger in turn inlines,
  69         -- so that if fromInteger is applied to an (S# i) the right thing happens
  70
  71     {-# INLINE fromInt #-}
  72     fromInt i           =  int2Float i
  73
  74 instance  Real Float  where
  75     toRational x        =  (m%1)*(b%1)^^n
  76                            where (m,n) = decodeFloat x
  77                                  b     = floatRadix  x
  78
  79 instance  Fractional Float  where
  80     (/) x y             =  divideFloat x y
  81     fromRational x      =  fromRat x
  82     recip x             =  1.0 / x
  83
  84 instance  Floating Float  where
  85     pi                  =  3.141592653589793238
  86     exp x               =  expFloat x
  87     log x               =  logFloat x
  88     sqrt x              =  sqrtFloat x
  89     sin x               =  sinFloat x
  90     cos x               =  cosFloat x
  91     tan x               =  tanFloat x
  92     asin x              =  asinFloat x
  93     acos x              =  acosFloat x
  94     atan x              =  atanFloat x
  95     sinh x              =  sinhFloat x
  96     cosh x              =  coshFloat x
  97     tanh x              =  tanhFloat x
  98     (**) x y            =  powerFloat x y
  99     logBase x y         =  log y / log x
 100
 101     asinh x = log (x + sqrt (1.0+x*x))
 102     acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
 103     atanh x = log ((x+1.0) / sqrt (1.0-x*x))
 104
 105 instance  RealFrac Float  where
 106
 107     {-# SPECIALIZE properFraction :: Float -> (Int, Float) #-}
 108     {-# SPECIALIZE truncate :: Float -> Int #-}
 109     {-# SPECIALIZE round    :: Float -> Int #-}
 110     {-# SPECIALIZE ceiling  :: Float -> Int #-}
 111     {-# SPECIALIZE floor    :: Float -> Int #-}
 112
 113     {-# SPECIALIZE properFraction :: Float -> (Integer, Float) #-}
 114     {-# SPECIALIZE truncate :: Float -> Integer #-}
 115     {-# SPECIALIZE round    :: Float -> Integer #-}
 116     {-# SPECIALIZE ceiling  :: Float -> Integer #-}
 117     {-# SPECIALIZE floor    :: Float -> Integer #-}
 118
 119     properFraction x
 120       = case (decodeFloat x)      of { (m,n) ->
 121         let  b = floatRadix x     in
 122         if n >= 0 then
 123             (fromInteger m * fromInteger b ^ n, 0.0)
 124         else
 125             case (quotRem m (b^(negate n))) of { (w,r) ->
 126             (fromInteger w, encodeFloat r n)
 127             }
 128         }
 129
 130     truncate x  = case properFraction x of
 131                      (n,_) -> n
 132
 133     round x     = case properFraction x of
 134                      (n,r) -> let
 135                                 m         = if r < 0.0 then n - 1 else n + 1
 136                                 half_down = abs r - 0.5
 137                               in
 138                               case (compare half_down 0.0) of
 139                                 LT -> n
 140                                 EQ -> if even n then n else m
 141                                 GT -> m
 142
 143     ceiling x   = case properFraction x of
 144                     (n,r) -> if r > 0.0 then n + 1 else n
 145
 146     floor x     = case properFraction x of
 147                     (n,r) -> if r < 0.0 then n - 1 else n
 148
 149 foreign import ccall "__encodeFloat" unsafe
 150         encodeFloat# :: Int# -> ByteArray# -> Int -> Float
 151 foreign import ccall "__int_encodeFloat" unsafe
 152         int_encodeFloat# :: Int# -> Int -> Float
 153
 154
 155 foreign import ccall "isFloatNaN" unsafe isFloatNaN :: Float -> Int
 156 foreign import ccall "isFloatInfinite" unsafe isFloatInfinite :: Float -> Int
 157 foreign import ccall "isFloatDenormalized" unsafe isFloatDenormalized :: Float -> Int
 158 foreign import ccall "isFloatNegativeZero" unsafe isFloatNegativeZero :: Float -> Int
 159
 160 instance  RealFloat Float  where
 161     floatRadix _        =  FLT_RADIX        -- from float.h
 162     floatDigits _       =  FLT_MANT_DIG     -- ditto
 163     floatRange _        =  (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
 164
 165     decodeFloat (F# f#)
 166       = case decodeFloat# f#    of
 167           (# exp#, s#, d# #) -> (J# s# d#, I# exp#)
 168
 169     encodeFloat (S# i) j     = int_encodeFloat# i j
 170     encodeFloat (J# s# d#) e = encodeFloat# s# d# e
 171
 172     exponent x          = case decodeFloat x of
 173                             (m,n) -> if m == 0 then 0 else n + floatDigits x
 174
 175     significand x       = case decodeFloat x of
 176                             (m,_) -> encodeFloat m (negate (floatDigits x))
 177
 178     scaleFloat k x      = case decodeFloat x of
 179                             (m,n) -> encodeFloat m (n+k)
 180     isNaN x          = 0 /= isFloatNaN x
 181     isInfinite x     = 0 /= isFloatInfinite x
 182     isDenormalized x = 0 /= isFloatDenormalized x
 183     isNegativeZero x = 0 /= isFloatNegativeZero x
 184     isIEEE _         = True
 185 \end{code}
 186
 187 %*********************************************************
 188 %*                                                      *
 189 \subsection{Type @Double@}
 190 %*                                                      *
 191 %*********************************************************
 192
 193 \begin{code}
 194 instance  Show Float  where
 195     showsPrec   x = showSigned showFloat x
 196     showList = showList__ (showsPrec 0)
 197
 198 instance Eq Double where
 199     (D# x) == (D# y) = x ==## y
 200
 201 instance Ord Double where
 202     (D# x) `compare` (D# y) | x <## y   = LT
 203                             | x ==## y  = EQ
 204                             | otherwise = GT
 205
 206     (D# x) <  (D# y) = x <##  y
 207     (D# x) <= (D# y) = x <=## y
 208     (D# x) >= (D# y) = x >=## y
 209     (D# x) >  (D# y) = x >##  y
 210
 211 instance  Num Double  where
 212     (+)         x y     =  plusDouble x y
 213     (-)         x y     =  minusDouble x y
 214     negate      x       =  negateDouble x
 215     (*)         x y     =  timesDouble x y
 216     abs x | x >= 0.0    =  x
 217           | otherwise   =  negateDouble x
 218     signum x | x == 0.0  = 0
 219              | x > 0.0   = 1
 220              | otherwise = negate 1
 221
 222     {-# INLINE fromInteger #-}
 223         -- See comments with Num Float
 224     fromInteger n       =  encodeFloat n 0
 225     fromInt (I# n#)     =  case (int2Double# n#) of { d# -> D# d# }
 226
 227 instance  Real Double  where
 228     toRational x        =  (m%1)*(b%1)^^n
 229                            where (m,n) = decodeFloat x
 230                                  b     = floatRadix  x
 231
 232 instance  Fractional Double  where
 233     (/) x y             =  divideDouble x y
 234     fromRational x      =  fromRat x
 235     recip x             =  1.0 / x
 236
 237 instance  Floating Double  where
 238     pi                  =  3.141592653589793238
 239     exp x               =  expDouble x
 240     log x               =  logDouble x
 241     sqrt x              =  sqrtDouble x
 242     sin  x              =  sinDouble x
 243     cos  x              =  cosDouble x
 244     tan  x              =  tanDouble x
 245     asin x              =  asinDouble x
 246     acos x              =  acosDouble x
 247     atan x              =  atanDouble x
 248     sinh x              =  sinhDouble x
 249     cosh x              =  coshDouble x
 250     tanh x              =  tanhDouble x
 251     (**) x y            =  powerDouble x y
 252     logBase x y         =  log y / log x
 253
 254     asinh x = log (x + sqrt (1.0+x*x))
 255     acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
 256     atanh x = log ((x+1.0) / sqrt (1.0-x*x))
 257
 258 instance  RealFrac Double  where
 259
 260     {-# SPECIALIZE properFraction :: Double -> (Int, Double) #-}
 261     {-# SPECIALIZE truncate :: Double -> Int #-}
 262     {-# SPECIALIZE round    :: Double -> Int #-}
 263     {-# SPECIALIZE ceiling  :: Double -> Int #-}
 264     {-# SPECIALIZE floor    :: Double -> Int #-}
 265
 266     {-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
 267     {-# SPECIALIZE truncate :: Double -> Integer #-}
 268     {-# SPECIALIZE round    :: Double -> Integer #-}
 269     {-# SPECIALIZE ceiling  :: Double -> Integer #-}
 270     {-# SPECIALIZE floor    :: Double -> Integer #-}
 271
 272     properFraction x
 273       = case (decodeFloat x)      of { (m,n) ->
 274         let  b = floatRadix x     in
 275         if n >= 0 then
 276             (fromInteger m * fromInteger b ^ n, 0.0)
 277         else
 278             case (quotRem m (b^(negate n))) of { (w,r) ->
 279             (fromInteger w, encodeFloat r n)
 280             }
 281         }
 282
 283     truncate x  = case properFraction x of
 284                      (n,_) -> n
 285
 286     round x     = case properFraction x of
 287                      (n,r) -> let
 288                                 m         = if r < 0.0 then n - 1 else n + 1
 289                                 half_down = abs r - 0.5
 290                               in
 291                               case (compare half_down 0.0) of
 292                                 LT -> n
 293                                 EQ -> if even n then n else m
 294                                 GT -> m
 295
 296     ceiling x   = case properFraction x of
 297                     (n,r) -> if r > 0.0 then n + 1 else n
 298
 299     floor x     = case properFraction x of
 300                     (n,r) -> if r < 0.0 then n - 1 else n
 301
 302 foreign import ccall "__encodeDouble" unsafe
 303         encodeDouble# :: Int# -> ByteArray# -> Int -> Double
 304 foreign import ccall "__int_encodeDouble" unsafe
 305         int_encodeDouble# :: Int# -> Int -> Double
 306
 307 foreign import ccall "isDoubleNaN" unsafe isDoubleNaN :: Double -> Int
 308 foreign import ccall "isDoubleInfinite" unsafe isDoubleInfinite :: Double -> Int
 309 foreign import ccall "isDoubleDenormalized" unsafe isDoubleDenormalized :: Double -> Int
 310 foreign import ccall "isDoubleNegativeZero" unsafe isDoubleNegativeZero :: Double -> Int
 311
 312 instance  RealFloat Double  where
 313     floatRadix _        =  FLT_RADIX        -- from float.h
 314     floatDigits _       =  DBL_MANT_DIG     -- ditto
 315     floatRange _        =  (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
 316
 317     decodeFloat (D# x#)
 318       = case decodeDouble# x#   of
 319           (# exp#, s#, d# #) -> (J# s# d#, I# exp#)
 320
 321     encodeFloat (S# i) j     = int_encodeDouble# i j
 322     encodeFloat (J# s# d#) e = encodeDouble# s# d# e
 323
 324     exponent x          = case decodeFloat x of
 325                             (m,n) -> if m == 0 then 0 else n + floatDigits x
 326
 327     significand x       = case decodeFloat x of
 328                             (m,_) -> encodeFloat m (negate (floatDigits x))
 329
 330     scaleFloat k x      = case decodeFloat x of
 331                             (m,n) -> encodeFloat m (n+k)
 332
 333     isNaN x             = 0 /= isDoubleNaN x
 334     isInfinite x        = 0 /= isDoubleInfinite x
 335     isDenormalized x    = 0 /= isDoubleDenormalized x
 336     isNegativeZero x    = 0 /= isDoubleNegativeZero x
 337     isIEEE _            = True
 338
 339 instance  Show Double  where
 340     showsPrec   x = showSigned showFloat x
 341     showList = showList__ (showsPrec 0)
 342 \end{code}
 343
 344 %*********************************************************
 345 %*                                                      *
 346 \subsection{Coercions}
 347 %*                                                      *
 348 %*********************************************************
 349
 350 \begin{code}
 351 {-# SPECIALIZE fromIntegral ::
 352     Int         -> Rational,
 353     Integer     -> Rational,
 354     Int         -> Int,
 355     Int         -> Integer,
 356     Int         -> Float,
 357     Int         -> Double,
 358     Integer     -> Int,
 359     Integer     -> Integer,
 360     Integer     -> Float,
 361     Integer     -> Double #-}
 362 fromIntegral    :: (Integral a, Num b) => a -> b
 363 fromIntegral    =  fromInteger . toInteger
 364
 365 {-# SPECIALIZE realToFrac ::
 366     Double      -> Rational,
 367     Rational    -> Double,
 368     Float       -> Rational,
 369     Rational    -> Float,
 370     Rational    -> Rational,
 371     Double      -> Double,
 372     Double      -> Float,
 373     Float       -> Float,
 374     Float       -> Double #-}
 375 realToFrac      :: (Real a, Fractional b) => a -> b
 376 realToFrac      =  fromRational . toRational
 377 \end{code}
 378
 379 %*********************************************************
 380 %*                                                      *
 381 \subsection{Common code for @Float@ and @Double@}
 382 %*                                                      *
 383 %*********************************************************
 384
 385 The @Enum@ instances for Floats and Doubles are slightly unusual.
 386 The @toEnum@ function truncates numbers to Int.  The definitions
 387 of @enumFrom@ and @enumFromThen@ allow floats to be used in arithmetic
 388 series: [0,0.1 .. 1.0].  However, roundoff errors make these somewhat
 389 dubious.  This example may have either 10 or 11 elements, depending on
 390 how 0.1 is represented.
 391
 392 NOTE: The instances for Float and Double do not make use of the default
 393 methods for @enumFromTo@ and @enumFromThenTo@, as these rely on there being
 394 a `non-lossy' conversion to and from Ints. Instead we make use of the
 395 1.2 default methods (back in the days when Enum had Ord as a superclass)
 396 for these (@numericEnumFromTo@ and @numericEnumFromThenTo@ below.)
 397
 398 \begin{code}
 399 instance  Enum Float  where
 400     succ x         = x + 1
 401     pred x         = x - 1
 402     toEnum         =  fromIntegral
 403     fromEnum       =  fromInteger . truncate   -- may overflow
 404     enumFrom       =  numericEnumFrom
 405     enumFromTo     =  numericEnumFromTo
 406     enumFromThen   =  numericEnumFromThen
 407     enumFromThenTo =  numericEnumFromThenTo
 408
 409 instance  Enum Double  where
 410     succ x         = x + 1
 411     pred x         = x - 1
 412     toEnum         =  fromIntegral
 413     fromEnum       =  fromInteger . truncate   -- may overflow
 414     enumFrom       =  numericEnumFrom
 415     enumFromTo     =  numericEnumFromTo
 416     enumFromThen   =  numericEnumFromThen
 417     enumFromThenTo =  numericEnumFromThenTo
 418
 419 numericEnumFrom         :: (Fractional a) => a -> [a]
 420 numericEnumFrom         =  iterate (+1)
 421
 422 numericEnumFromThen     :: (Fractional a) => a -> a -> [a]
 423 numericEnumFromThen n m =  iterate (+(m-n)) n
 424
 425 numericEnumFromTo       :: (Ord a, Fractional a) => a -> a -> [a]
 426 numericEnumFromTo n m   = takeWhile (<= m + 1/2) (numericEnumFrom n)
 427
 428 numericEnumFromThenTo   :: (Ord a, Fractional a) => a -> a -> a -> [a]
 429 numericEnumFromThenTo e1 e2 e3 = takeWhile pred (numericEnumFromThen e1 e2)
 430                                 where
 431                                  mid = (e2 - e1) / 2
 432                                  pred | e2 > e1   = (<= e3 + mid)
 433                                       | otherwise = (>= e3 + mid)
 434
 435 \end{code}
 436
 437 @approxRational@, applied to two real fractional numbers x and epsilon,
 438 returns the simplest rational number within epsilon of x.  A rational
 439 number n%d in reduced form is said to be simpler than another n'%d' if
 440 abs n <= abs n' && d <= d'.  Any real interval contains a unique
 441 simplest rational; here, for simplicity, we assume a closed rational
 442 interval.  If such an interval includes at least one whole number, then
 443 the simplest rational is the absolutely least whole number.  Otherwise,
 444 the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
 445 and abs r' < d', and the simplest rational is q%1 + the reciprocal of
 446 the simplest rational between d'%r' and d%r.
 447
 448 \begin{code}
 449 approxRational          :: (RealFrac a) => a -> a -> Rational
 450 approxRational rat eps  =  simplest (rat-eps) (rat+eps)
 451         where simplest x y | y < x      =  simplest y x
 452                            | x == y     =  xr
 453                            | x > 0      =  simplest' n d n' d'
 454                            | y < 0      =  - simplest' (-n') d' (-n) d
 455                            | otherwise  =  0 :% 1
 456                                         where xr  = toRational x
 457                                               n   = numerator xr
 458                                               d   = denominator xr
 459                                               nd' = toRational y
 460                                               n'  = numerator nd'
 461                                               d'  = denominator nd'
 462
 463               simplest' n d n' d'       -- assumes 0 < n%d < n'%d'
 464                         | r == 0     =  q :% 1
 465                         | q /= q'    =  (q+1) :% 1
 466                         | otherwise  =  (q*n''+d'') :% n''
 467                                      where (q,r)      =  quotRem n d
 468                                            (q',r')    =  quotRem n' d'
 469                                            nd''       =  simplest' d' r' d r
 470                                            n''        =  numerator nd''
 471                                            d''        =  denominator nd''
 472 \end{code}
 473
 474
 475 \begin{code}
 476 instance  (Integral a)  => Ord (Ratio a)  where
 477     (x:%y) <= (x':%y')  =  x * y' <= x' * y
 478     (x:%y) <  (x':%y')  =  x * y' <  x' * y
 479
 480 instance  (Integral a)  => Num (Ratio a)  where
 481     (x:%y) + (x':%y')   =  reduce (x*y' + x'*y) (y*y')
 482     (x:%y) - (x':%y')   =  reduce (x*y' - x'*y) (y*y')
 483     (x:%y) * (x':%y')   =  reduce (x * x') (y * y')
 484     negate (x:%y)       =  (-x) :% y
 485     abs (x:%y)          =  abs x :% y
 486     signum (x:%_)       =  signum x :% 1
 487     fromInteger x       =  fromInteger x :% 1
 488
 489 instance  (Integral a)  => Real (Ratio a)  where
 490     toRational (x:%y)   =  toInteger x :% toInteger y
 491
 492 instance  (Integral a)  => Fractional (Ratio a)  where
 493     (x:%y) / (x':%y')   =  (x*y') % (y*x')
 494     recip (x:%y)        =  if x < 0 then (-y) :% (-x) else y :% x
 495     fromRational (x:%y) =  fromInteger x :% fromInteger y
 496
 497 instance  (Integral a)  => RealFrac (Ratio a)  where
 498     properFraction (x:%y) = (fromIntegral q, r:%y)
 499                             where (q,r) = quotRem x y
 500
 501 instance  (Integral a)  => Enum (Ratio a)  where
 502     succ x              =  x + 1
 503     pred x              =  x - 1
 504
 505     toEnum n            =  fromIntegral n :% 1
 506     fromEnum            =  fromInteger . truncate
 507
 508     enumFrom            =  bounded_iterator True (1)
 509     enumFromThen n m    =  bounded_iterator (diff >= 0) diff n
 510                           where diff = m - n
 511
 512
 513 bounded_iterator :: (Ord a, Num a) => Bool -> a -> a -> [a]
 514 bounded_iterator inc step v
 515    | inc      && v > new_v = [v]  -- oflow
 516    | not inc  && v < new_v = [v]  -- uflow
 517    | otherwise             = v : bounded_iterator inc step new_v
 518   where
 519    new_v = v + step
 520
 521 ratio_prec :: Int
 522 ratio_prec = 7
 523
 524 instance  (Integral a)  => Show (Ratio a)  where
 525     showsPrec p (x:%y)  =  showParen (p > ratio_prec)
 526                                (shows x . showString " % " . shows y)
 527 \end{code}
 528
 529 @showRational@ converts a Rational to a string that looks like a
 530 floating point number, but without converting to any floating type
 531 (because of the possible overflow).
 532
 533 From/by Lennart, 94/09/26
 534
 535 \begin{code}
 536 showRational :: Int -> Rational -> String
 537 showRational n r =
 538     if r == 0 then
 539         "0.0"
 540     else
 541         let (r', e) = normalize r
 542         in  prR n r' e
 543
 544 startExpExp :: Int
 545 startExpExp = 4
 546
 547 -- make sure 1 <= r < 10
 548 normalize :: Rational -> (Rational, Int)
 549 normalize r = if r < 1 then
 550                   case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
 551               else
 552                   norm startExpExp r 0
 553         where norm :: Int -> Rational -> Int -> (Rational, Int)
 554               -- Invariant: x*10^e == original r
 555               norm 0  x e = (x, e)
 556               norm ee x e =
 557                 let n = 10^ee
 558                     tn = 10^n
 559                 in  if x >= tn then norm ee (x/tn) (e+n) else norm (ee-1) x e
 560
 561 prR :: Int -> Rational -> Int -> String
 562 prR n r e  | r <  1  = prR n (r*10) (e-1)               -- final adjustment
 563 prR n r e  | r >= 10 = prR n (r/10) (e+1)
 564 prR n r e0
 565   | e > 0 && e < 8   = takeN e s ('.' : drop0 (drop e s) [])
 566   | e <= 0 && e > -3 = '0': '.' : takeN (-e) (repeat '0') (drop0 s [])
 567   | otherwise        =  h : '.' : drop0 t ('e':show e0)
 568    where
 569         s@(h:t) = show ((round (r * 10^n))::Integer)
 570         e       = e0+1
 571
 572 #ifdef USE_REPORT_PRELUDE
 573         takeN n ls rs = take n ls ++ rs
 574 #else
 575         takeN (I# n#) ls rs = takeUInt_append n# ls rs
 576 #endif
 577
 578 drop0 :: String -> String -> String
 579 drop0     [] rs = rs
 580 drop0 (c:cs) rs = c : fromMaybe rs (dropTrailing0s cs) --WAS (yuck): reverse (dropWhile (=='0') (reverse cs))
 581   where
 582    dropTrailing0s []       = Nothing
 583    dropTrailing0s ('0':xs) =
 584      case dropTrailing0s xs of
 585        Nothing -> Nothing
 586        Just ls -> Just ('0':ls)
 587    dropTrailing0s (x:xs) =
 588      case dropTrailing0s xs of
 589       Nothing -> Just [x]
 590       Just ls -> Just (x:ls)
 591
 592 \end{code}
 593
 594 [In response to a request for documentation of how fromRational works,
 595 Joe Fasel writes:] A quite reasonable request!  This code was added to
 596 the Prelude just before the 1.2 release, when Lennart, working with an
 597 early version of hbi, noticed that (read . show) was not the identity
 598 for floating-point numbers.  (There was a one-bit error about half the
 599 time.)  The original version of the conversion function was in fact
 600 simply a floating-point divide, as you suggest above. The new version
 601 is, I grant you, somewhat denser.
 602
 603 Unfortunately, Joe's code doesn't work!  Here's an example:
 604
 605 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
 606
 607 This program prints
 608         0.0000000000000000
 609 instead of
 610         1.8217369128763981e-300
 611
 612 Lennart's code follows, and it works...
 613
 614 \begin{pseudocode}
 615 fromRat :: (RealFloat a) => Rational -> a
 616 fromRat x = x'
 617         where x' = f e
 618
 619 --              If the exponent of the nearest floating-point number to x
 620 --              is e, then the significand is the integer nearest xb^(-e),
 621 --              where b is the floating-point radix.  We start with a good
 622 --              guess for e, and if it is correct, the exponent of the
 623 --              floating-point number we construct will again be e.  If
 624 --              not, one more iteration is needed.
 625
 626               f e   = if e' == e then y else f e'
 627                       where y      = encodeFloat (round (x * (1 % b)^^e)) e
 628                             (_,e') = decodeFloat y
 629               b     = floatRadix x'
 630
 631 --              We obtain a trial exponent by doing a floating-point
 632 --              division of x's numerator by its denominator.  The
 633 --              result of this division may not itself be the ultimate
 634 --              result, because of an accumulation of three rounding
 635 --              errors.
 636
 637               (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
 638                                         / fromInteger (denominator x))
 639 \end{pseudocode}
 640
 641 Now, here's Lennart's code.
 642
 643 \begin{code}
 644 {-# SPECIALISE fromRat ::
 645         Rational -> Double,
 646         Rational -> Float #-}
 647 fromRat :: (RealFloat a) => Rational -> a
 648 fromRat x
 649   | x == 0    =  encodeFloat 0 0                -- Handle exceptional cases
 650   | x <  0    =  - fromRat' (-x)                -- first.
 651   | otherwise =  fromRat' x
 652
 653 -- Conversion process:
 654 -- Scale the rational number by the RealFloat base until
 655 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
 656 -- Then round the rational to an Integer and encode it with the exponent
 657 -- that we got from the scaling.
 658 -- To speed up the scaling process we compute the log2 of the number to get
 659 -- a first guess of the exponent.
 660
 661 fromRat' :: (RealFloat a) => Rational -> a
 662 fromRat' x = r
 663   where b = floatRadix r
 664         p = floatDigits r
 665         (minExp0, _) = floatRange r
 666         minExp = minExp0 - p            -- the real minimum exponent
 667         xMin   = toRational (expt b (p-1))
 668         xMax   = toRational (expt b p)
 669         p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
 670         f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
 671         (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
 672         r = encodeFloat (round x') p'
 673
 674 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
 675 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
 676 scaleRat b minExp xMin xMax p x
 677  | p <= minExp = (x, p)
 678  | x >= xMax   = scaleRat b minExp xMin xMax (p+1) (x/b)
 679  | x < xMin    = scaleRat b minExp xMin xMax (p-1) (x*b)
 680  | otherwise   = (x, p)
 681
 682 -- Exponentiation with a cache for the most common numbers.
 683 minExpt, maxExpt :: Int
 684 minExpt = 0
 685 maxExpt = 1100
 686
 687 expt :: Integer -> Int -> Integer
 688 expt base n =
 689     if base == 2 && n >= minExpt && n <= maxExpt then
 690         expts!n
 691     else
 692         base^n
 693
 694 expts :: Array Int Integer
 695 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
 696
 697 -- Compute the (floor of the) log of i in base b.
 698 -- Simplest way would be just divide i by b until it's smaller then b, but that would
 699 -- be very slow!  We are just slightly more clever.
 700 integerLogBase :: Integer -> Integer -> Int
 701 integerLogBase b i
 702    | i < b     = 0
 703    | otherwise = doDiv (i `div` (b^l)) l
 704        where
 705         -- Try squaring the base first to cut down the number of divisions.
 706          l = 2 * integerLogBase (b*b) i
 707
 708          doDiv :: Integer -> Int -> Int
 709          doDiv x y
 710             | x < b     = y
 711             | otherwise = doDiv (x `div` b) (y+1)
 712
 713 \end{code}
 714
 715 %*********************************************************
 716 %*                                                      *
 717 \subsection{Printing out numbers}
 718 %*                                                      *
 719 %*********************************************************
 720
 721 \begin{code}
 722 --Exported from std library Numeric, defined here to
 723 --avoid mut. rec. between PrelNum and Numeric.
 724 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
 725 showSigned showPos p x
 726    | x < 0     = showParen (p > 6) (showChar '-' . showPos (-x))
 727    | otherwise = showPos x
 728
 729 showFloat :: (RealFloat a) => a -> ShowS
 730 showFloat x  =  showString (formatRealFloat FFGeneric Nothing x)
 731
 732 -- These are the format types.  This type is not exported.
 733
 734 data FFFormat = FFExponent | FFFixed | FFGeneric
 735
 736 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
 737 formatRealFloat fmt decs x
 738    | isNaN x                   = "NaN"
 739    | isInfinite x              = if x < 0 then "-Infinity" else "Infinity"
 740    | x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
 741    | otherwise                 = doFmt fmt (floatToDigits (toInteger base) x)
 742  where
 743   base = 10
 744
 745   doFmt format (is, e) =
 746     let ds = map intToDigit is in
 747     case format of
 748      FFGeneric ->
 749       doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
 750             (is,e)
 751      FFExponent ->
 752       case decs of
 753        Nothing ->
 754         let e' = if e==0 then 0 else e-1 in
 755         (case ds of
 756           [d]     -> d : ".0e"
 757           (d:ds') -> d : '.' : ds' ++ "e") ++ show e'
 758        Just dec ->
 759         let dec' = max dec 1 in
 760         case is of
 761          [0] -> '0' :'.' : take dec' (repeat '0') ++ "e0"
 762          _ ->
 763           let
 764            (ei,is') = roundTo base (dec'+1) is
 765            (d:ds') = map intToDigit (if ei > 0 then init is' else is')
 766           in
 767           d:'.':ds' ++ 'e':show (e-1+ei)
 768      FFFixed ->
 769       let
 770        mk0 ls = case ls of { "" -> "0" ; _ -> ls}
 771       in
 772       case decs of
 773        Nothing ->
 774          let
 775           f 0 s    rs  = mk0 (reverse s) ++ '.':mk0 rs
 776           f n s    ""  = f (n-1) ('0':s) ""
 777           f n s (r:rs) = f (n-1) (r:s) rs
 778          in
 779          f e "" ds
 780        Just dec ->
 781         let dec' = max dec 0 in
 782         if e >= 0 then
 783          let
 784           (ei,is') = roundTo base (dec' + e) is
 785           (ls,rs)  = splitAt (e+ei) (map intToDigit is')
 786          in
 787          mk0 ls ++ (if null rs then "" else '.':rs)
 788         else
 789          let
 790           (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
 791           d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
 792          in
 793          d : '.' : ds'
 794
 795
 796 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
 797 roundTo base d is =
 798   case f d is of
 799     x@(0,_) -> x
 800     (1,xs)  -> (1, 1:xs)
 801  where
 802   b2 = base `div` 2
 803
 804   f n []     = (0, replicate n 0)
 805   f 0 (x:_)  = (if x >= b2 then 1 else 0, [])
 806   f n (i:xs)
 807      | i' == base = (1,0:ds)
 808      | otherwise  = (0,i':ds)
 809       where
 810        (c,ds) = f (n-1) xs
 811        i'     = c + i
 812
 813 --
 814 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
 815 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
 816 -- This version uses a much slower logarithm estimator. It should be improved.
 817
 818 -- This function returns a list of digits (Ints in [0..base-1]) and an
 819 -- exponent.
 820
 821 floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
 822 floatToDigits _ 0 = ([0], 0)
 823 floatToDigits base x =
 824  let
 825   (f0, e0) = decodeFloat x
 826   (minExp0, _) = floatRange x
 827   p = floatDigits x
 828   b = floatRadix x
 829   minExp = minExp0 - p -- the real minimum exponent
 830   -- Haskell requires that f be adjusted so denormalized numbers
 831   -- will have an impossibly low exponent.  Adjust for this.
 832   (f, e) =
 833    let n = minExp - e0 in
 834    if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
 835   (r, s, mUp, mDn) =
 836    if e >= 0 then
 837     let be = b^ e in
 838     if f == b^(p-1) then
 839       (f*be*b*2, 2*b, be*b, b)
 840     else
 841       (f*be*2, 2, be, be)
 842    else
 843     if e > minExp && f == b^(p-1) then
 844       (f*b*2, b^(-e+1)*2, b, 1)
 845     else
 846       (f*2, b^(-e)*2, 1, 1)
 847   k =
 848    let
 849     k0 =
 850      if b == 2 && base == 10 then
 851         -- logBase 10 2 is slightly bigger than 3/10 so
 852         -- the following will err on the low side.  Ignoring
 853         -- the fraction will make it err even more.
 854         -- Haskell promises that p-1 <= logBase b f < p.
 855         (p - 1 + e0) * 3 `div` 10
 856      else
 857         ceiling ((log (fromInteger (f+1)) +
 858                  fromInt e * log (fromInteger b)) /
 859                    log (fromInteger base))
 860 --WAS:            fromInt e * log (fromInteger b))
 861
 862     fixup n =
 863       if n >= 0 then
 864         if r + mUp <= expt base n * s then n else fixup (n+1)
 865       else
 866         if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
 867    in
 868    fixup k0
 869
 870   gen ds rn sN mUpN mDnN =
 871    let
 872     (dn, rn') = (rn * base) `divMod` sN
 873     mUpN' = mUpN * base
 874     mDnN' = mDnN * base
 875    in
 876    case (rn' < mDnN', rn' + mUpN' > sN) of
 877     (True,  False) -> dn : ds
 878     (False, True)  -> dn+1 : ds
 879     (True,  True)  -> if rn' * 2 < sN then dn : ds else dn+1 : ds
 880     (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
 881
 882   rds =
 883    if k >= 0 then
 884       gen [] r (s * expt base k) mUp mDn
 885    else
 886      let bk = expt base (-k) in
 887      gen [] (r * bk) s (mUp * bk) (mDn * bk)
 888  in
 889  (map toInt (reverse rds), k)
 890
 891 \end{code}
 892
 893 %*********************************************************
 894 %*                                                      *
 895 \subsection{Numeric primops}
 896 %*                                                      *
 897 %*********************************************************
 898
 899 Definitions of the boxed PrimOps; these will be
 900 used in the case of partial applications, etc.
 901
 902 \begin{code}
 903 plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
 904 plusFloat   (F# x) (F# y) = F# (plusFloat# x y)
 905 minusFloat  (F# x) (F# y) = F# (minusFloat# x y)
 906 timesFloat  (F# x) (F# y) = F# (timesFloat# x y)
 907 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
 908
 909 negateFloat :: Float -> Float
 910 negateFloat (F# x)        = F# (negateFloat# x)
 911
 912 gtFloat, geFloat, eqFloat, neFloat, ltFloat, leFloat :: Float -> Float -> Bool
 913 gtFloat     (F# x) (F# y) = gtFloat# x y
 914 geFloat     (F# x) (F# y) = geFloat# x y
 915 eqFloat     (F# x) (F# y) = eqFloat# x y
 916 neFloat     (F# x) (F# y) = neFloat# x y
 917 ltFloat     (F# x) (F# y) = ltFloat# x y
 918 leFloat     (F# x) (F# y) = leFloat# x y
 919
 920 float2Int :: Float -> Int
 921 float2Int   (F# x) = I# (float2Int# x)
 922
 923 int2Float :: Int -> Float
 924 int2Float   (I# x) = F# (int2Float# x)
 925
 926 expFloat, logFloat, sqrtFloat :: Float -> Float
 927 sinFloat, cosFloat, tanFloat  :: Float -> Float
 928 asinFloat, acosFloat, atanFloat  :: Float -> Float
 929 sinhFloat, coshFloat, tanhFloat  :: Float -> Float
 930 expFloat    (F# x) = F# (expFloat# x)
 931 logFloat    (F# x) = F# (logFloat# x)
 932 sqrtFloat   (F# x) = F# (sqrtFloat# x)
 933 sinFloat    (F# x) = F# (sinFloat# x)
 934 cosFloat    (F# x) = F# (cosFloat# x)
 935 tanFloat    (F# x) = F# (tanFloat# x)
 936 asinFloat   (F# x) = F# (asinFloat# x)
 937 acosFloat   (F# x) = F# (acosFloat# x)
 938 atanFloat   (F# x) = F# (atanFloat# x)
 939 sinhFloat   (F# x) = F# (sinhFloat# x)
 940 coshFloat   (F# x) = F# (coshFloat# x)
 941 tanhFloat   (F# x) = F# (tanhFloat# x)
 942
 943 powerFloat :: Float -> Float -> Float
 944 powerFloat  (F# x) (F# y) = F# (powerFloat# x y)
 945
 946 -- definitions of the boxed PrimOps; these will be
 947 -- used in the case of partial applications, etc.
 948
 949 plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
 950 plusDouble   (D# x) (D# y) = D# (x +## y)
 951 minusDouble  (D# x) (D# y) = D# (x -## y)
 952 timesDouble  (D# x) (D# y) = D# (x *## y)
 953 divideDouble (D# x) (D# y) = D# (x /## y)
 954
 955 negateDouble :: Double -> Double
 956 negateDouble (D# x)        = D# (negateDouble# x)
 957
 958 gtDouble, geDouble, eqDouble, neDouble, leDouble, ltDouble :: Double -> Double -> Bool
 959 gtDouble    (D# x) (D# y) = x >## y
 960 geDouble    (D# x) (D# y) = x >=## y
 961 eqDouble    (D# x) (D# y) = x ==## y
 962 neDouble    (D# x) (D# y) = x /=## y
 963 ltDouble    (D# x) (D# y) = x <## y
 964 leDouble    (D# x) (D# y) = x <=## y
 965
 966 double2Int :: Double -> Int
 967 double2Int   (D# x) = I# (double2Int#   x)
 968
 969 int2Double :: Int -> Double
 970 int2Double   (I# x) = D# (int2Double#   x)
 971
 972 double2Float :: Double -> Float
 973 double2Float (D# x) = F# (double2Float# x)
 974 float2Double :: Float -> Double
 975 float2Double (F# x) = D# (float2Double# x)
 976
 977 expDouble, logDouble, sqrtDouble :: Double -> Double
 978 sinDouble, cosDouble, tanDouble  :: Double -> Double
 979 asinDouble, acosDouble, atanDouble  :: Double -> Double
 980 sinhDouble, coshDouble, tanhDouble  :: Double -> Double
 981 expDouble    (D# x) = D# (expDouble# x)
 982 logDouble    (D# x) = D# (logDouble# x)
 983 sqrtDouble   (D# x) = D# (sqrtDouble# x)
 984 sinDouble    (D# x) = D# (sinDouble# x)
 985 cosDouble    (D# x) = D# (cosDouble# x)
 986 tanDouble    (D# x) = D# (tanDouble# x)
 987 asinDouble   (D# x) = D# (asinDouble# x)
 988 acosDouble   (D# x) = D# (acosDouble# x)
 989 atanDouble   (D# x) = D# (atanDouble# x)
 990 sinhDouble   (D# x) = D# (sinhDouble# x)
 991 coshDouble   (D# x) = D# (coshDouble# x)
 992 tanhDouble   (D# x) = D# (tanhDouble# x)
 993
 994 powerDouble :: Double -> Double -> Double
 995 powerDouble  (D# x) (D# y) = D# (x **## y)
 996 \end{code}