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         takeN (I# n#) ls rs = takeUInt_append n# ls rs
 573
 574 drop0 :: String -> String -> String
 575 drop0     [] rs = rs
 576 drop0 (c:cs) rs = c : fromMaybe rs (dropTrailing0s cs) --WAS (yuck): reverse (dropWhile (=='0') (reverse cs))
 577   where
 578    dropTrailing0s []       = Nothing
 579    dropTrailing0s ('0':xs) =
 580      case dropTrailing0s xs of
 581        Nothing -> Nothing
 582        Just ls -> Just ('0':ls)
 583    dropTrailing0s (x:xs) =
 584      case dropTrailing0s xs of
 585       Nothing -> Just [x]
 586       Just ls -> Just (x:ls)
 587
 588 \end{code}
 589
 590 [In response to a request for documentation of how fromRational works,
 591 Joe Fasel writes:] A quite reasonable request!  This code was added to
 592 the Prelude just before the 1.2 release, when Lennart, working with an
 593 early version of hbi, noticed that (read . show) was not the identity
 594 for floating-point numbers.  (There was a one-bit error about half the
 595 time.)  The original version of the conversion function was in fact
 596 simply a floating-point divide, as you suggest above. The new version
 597 is, I grant you, somewhat denser.
 598
 599 Unfortunately, Joe's code doesn't work!  Here's an example:
 600
 601 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
 602
 603 This program prints
 604         0.0000000000000000
 605 instead of
 606         1.8217369128763981e-300
 607
 608 Lennart's code follows, and it works...
 609
 610 \begin{pseudocode}
 611 fromRat :: (RealFloat a) => Rational -> a
 612 fromRat x = x'
 613         where x' = f e
 614
 615 --              If the exponent of the nearest floating-point number to x
 616 --              is e, then the significand is the integer nearest xb^(-e),
 617 --              where b is the floating-point radix.  We start with a good
 618 --              guess for e, and if it is correct, the exponent of the
 619 --              floating-point number we construct will again be e.  If
 620 --              not, one more iteration is needed.
 621
 622               f e   = if e' == e then y else f e'
 623                       where y      = encodeFloat (round (x * (1 % b)^^e)) e
 624                             (_,e') = decodeFloat y
 625               b     = floatRadix x'
 626
 627 --              We obtain a trial exponent by doing a floating-point
 628 --              division of x's numerator by its denominator.  The
 629 --              result of this division may not itself be the ultimate
 630 --              result, because of an accumulation of three rounding
 631 --              errors.
 632
 633               (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
 634                                         / fromInteger (denominator x))
 635 \end{pseudocode}
 636
 637 Now, here's Lennart's code.
 638
 639 \begin{code}
 640 {-# SPECIALISE fromRat ::
 641         Rational -> Double,
 642         Rational -> Float #-}
 643 fromRat :: (RealFloat a) => Rational -> a
 644 fromRat x
 645   | x == 0    =  encodeFloat 0 0                -- Handle exceptional cases
 646   | x <  0    =  - fromRat' (-x)                -- first.
 647   | otherwise =  fromRat' x
 648
 649 -- Conversion process:
 650 -- Scale the rational number by the RealFloat base until
 651 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
 652 -- Then round the rational to an Integer and encode it with the exponent
 653 -- that we got from the scaling.
 654 -- To speed up the scaling process we compute the log2 of the number to get
 655 -- a first guess of the exponent.
 656
 657 fromRat' :: (RealFloat a) => Rational -> a
 658 fromRat' x = r
 659   where b = floatRadix r
 660         p = floatDigits r
 661         (minExp0, _) = floatRange r
 662         minExp = minExp0 - p            -- the real minimum exponent
 663         xMin   = toRational (expt b (p-1))
 664         xMax   = toRational (expt b p)
 665         p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
 666         f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
 667         (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
 668         r = encodeFloat (round x') p'
 669
 670 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
 671 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
 672 scaleRat b minExp xMin xMax p x
 673  | p <= minExp = (x, p)
 674  | x >= xMax   = scaleRat b minExp xMin xMax (p+1) (x/b)
 675  | x < xMin    = scaleRat b minExp xMin xMax (p-1) (x*b)
 676  | otherwise   = (x, p)
 677
 678 -- Exponentiation with a cache for the most common numbers.
 679 minExpt, maxExpt :: Int
 680 minExpt = 0
 681 maxExpt = 1100
 682
 683 expt :: Integer -> Int -> Integer
 684 expt base n =
 685     if base == 2 && n >= minExpt && n <= maxExpt then
 686         expts!n
 687     else
 688         base^n
 689
 690 expts :: Array Int Integer
 691 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
 692
 693 -- Compute the (floor of the) log of i in base b.
 694 -- Simplest way would be just divide i by b until it's smaller then b, but that would
 695 -- be very slow!  We are just slightly more clever.
 696 integerLogBase :: Integer -> Integer -> Int
 697 integerLogBase b i
 698    | i < b     = 0
 699    | otherwise = doDiv (i `div` (b^l)) l
 700        where
 701         -- Try squaring the base first to cut down the number of divisions.
 702          l = 2 * integerLogBase (b*b) i
 703
 704          doDiv :: Integer -> Int -> Int
 705          doDiv x y
 706             | x < b     = y
 707             | otherwise = doDiv (x `div` b) (y+1)
 708
 709 \end{code}
 710
 711 %*********************************************************
 712 %*                                                      *
 713 \subsection{Printing out numbers}
 714 %*                                                      *
 715 %*********************************************************
 716
 717 \begin{code}
 718 --Exported from std library Numeric, defined here to
 719 --avoid mut. rec. between PrelNum and Numeric.
 720 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
 721 showSigned showPos p x
 722    | x < 0     = showParen (p > 6) (showChar '-' . showPos (-x))
 723    | otherwise = showPos x
 724
 725 showFloat :: (RealFloat a) => a -> ShowS
 726 showFloat x  =  showString (formatRealFloat FFGeneric Nothing x)
 727
 728 -- These are the format types.  This type is not exported.
 729
 730 data FFFormat = FFExponent | FFFixed | FFGeneric
 731
 732 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
 733 formatRealFloat fmt decs x
 734    | isNaN x                   = "NaN"
 735    | isInfinite x              = if x < 0 then "-Infinity" else "Infinity"
 736    | x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
 737    | otherwise                 = doFmt fmt (floatToDigits (toInteger base) x)
 738  where
 739   base = 10
 740
 741   doFmt format (is, e) =
 742     let ds = map intToDigit is in
 743     case format of
 744      FFGeneric ->
 745       doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
 746             (is,e)
 747      FFExponent ->
 748       case decs of
 749        Nothing ->
 750         let e' = if e==0 then 0 else e-1 in
 751         (case ds of
 752           [d]     -> d : ".0e"
 753           (d:ds') -> d : '.' : ds' ++ "e") ++ show e'
 754        Just dec ->
 755         let dec' = max dec 1 in
 756         case is of
 757          [0] -> '0' :'.' : take dec' (repeat '0') ++ "e0"
 758          _ ->
 759           let
 760            (ei,is') = roundTo base (dec'+1) is
 761            (d:ds') = map intToDigit (if ei > 0 then init is' else is')
 762           in
 763           d:'.':ds' ++ 'e':show (e-1+ei)
 764      FFFixed ->
 765       let
 766        mk0 ls = case ls of { "" -> "0" ; _ -> ls}
 767       in
 768       case decs of
 769        Nothing ->
 770          let
 771           f 0 s    rs  = mk0 (reverse s) ++ '.':mk0 rs
 772           f n s    ""  = f (n-1) ('0':s) ""
 773           f n s (r:rs) = f (n-1) (r:s) rs
 774          in
 775          f e "" ds
 776        Just dec ->
 777         let dec' = max dec 0 in
 778         if e >= 0 then
 779          let
 780           (ei,is') = roundTo base (dec' + e) is
 781           (ls,rs)  = splitAt (e+ei) (map intToDigit is')
 782          in
 783          mk0 ls ++ (if null rs then "" else '.':rs)
 784         else
 785          let
 786           (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
 787           d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
 788          in
 789          d : '.' : ds'
 790
 791
 792 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
 793 roundTo base d is =
 794   case f d is of
 795     x@(0,_) -> x
 796     (1,xs)  -> (1, 1:xs)
 797  where
 798   b2 = base `div` 2
 799
 800   f n []     = (0, replicate n 0)
 801   f 0 (x:_)  = (if x >= b2 then 1 else 0, [])
 802   f n (i:xs)
 803      | i' == base = (1,0:ds)
 804      | otherwise  = (0,i':ds)
 805       where
 806        (c,ds) = f (n-1) xs
 807        i'     = c + i
 808
 809 --
 810 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
 811 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
 812 -- This version uses a much slower logarithm estimator. It should be improved.
 813
 814 -- This function returns a list of digits (Ints in [0..base-1]) and an
 815 -- exponent.
 816
 817 floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
 818 floatToDigits _ 0 = ([0], 0)
 819 floatToDigits base x =
 820  let
 821   (f0, e0) = decodeFloat x
 822   (minExp0, _) = floatRange x
 823   p = floatDigits x
 824   b = floatRadix x
 825   minExp = minExp0 - p -- the real minimum exponent
 826   -- Haskell requires that f be adjusted so denormalized numbers
 827   -- will have an impossibly low exponent.  Adjust for this.
 828   (f, e) =
 829    let n = minExp - e0 in
 830    if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
 831   (r, s, mUp, mDn) =
 832    if e >= 0 then
 833     let be = b^ e in
 834     if f == b^(p-1) then
 835       (f*be*b*2, 2*b, be*b, b)
 836     else
 837       (f*be*2, 2, be, be)
 838    else
 839     if e > minExp && f == b^(p-1) then
 840       (f*b*2, b^(-e+1)*2, b, 1)
 841     else
 842       (f*2, b^(-e)*2, 1, 1)
 843   k =
 844    let
 845     k0 =
 846      if b == 2 && base == 10 then
 847         -- logBase 10 2 is slightly bigger than 3/10 so
 848         -- the following will err on the low side.  Ignoring
 849         -- the fraction will make it err even more.
 850         -- Haskell promises that p-1 <= logBase b f < p.
 851         (p - 1 + e0) * 3 `div` 10
 852      else
 853         ceiling ((log (fromInteger (f+1)) +
 854                  fromInt e * log (fromInteger b)) /
 855                    log (fromInteger base))
 856 --WAS:            fromInt e * log (fromInteger b))
 857
 858     fixup n =
 859       if n >= 0 then
 860         if r + mUp <= expt base n * s then n else fixup (n+1)
 861       else
 862         if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
 863    in
 864    fixup k0
 865
 866   gen ds rn sN mUpN mDnN =
 867    let
 868     (dn, rn') = (rn * base) `divMod` sN
 869     mUpN' = mUpN * base
 870     mDnN' = mDnN * base
 871    in
 872    case (rn' < mDnN', rn' + mUpN' > sN) of
 873     (True,  False) -> dn : ds
 874     (False, True)  -> dn+1 : ds
 875     (True,  True)  -> if rn' * 2 < sN then dn : ds else dn+1 : ds
 876     (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
 877
 878   rds =
 879    if k >= 0 then
 880       gen [] r (s * expt base k) mUp mDn
 881    else
 882      let bk = expt base (-k) in
 883      gen [] (r * bk) s (mUp * bk) (mDn * bk)
 884  in
 885  (map toInt (reverse rds), k)
 886
 887 \end{code}
 888
 889 %*********************************************************
 890 %*                                                      *
 891 \subsection{Numeric primops}
 892 %*                                                      *
 893 %*********************************************************
 894
 895 Definitions of the boxed PrimOps; these will be
 896 used in the case of partial applications, etc.
 897
 898 \begin{code}
 899 plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
 900 plusFloat   (F# x) (F# y) = F# (plusFloat# x y)
 901 minusFloat  (F# x) (F# y) = F# (minusFloat# x y)
 902 timesFloat  (F# x) (F# y) = F# (timesFloat# x y)
 903 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
 904
 905 negateFloat :: Float -> Float
 906 negateFloat (F# x)        = F# (negateFloat# x)
 907
 908 gtFloat, geFloat, eqFloat, neFloat, ltFloat, leFloat :: Float -> Float -> Bool
 909 gtFloat     (F# x) (F# y) = gtFloat# x y
 910 geFloat     (F# x) (F# y) = geFloat# x y
 911 eqFloat     (F# x) (F# y) = eqFloat# x y
 912 neFloat     (F# x) (F# y) = neFloat# x y
 913 ltFloat     (F# x) (F# y) = ltFloat# x y
 914 leFloat     (F# x) (F# y) = leFloat# x y
 915
 916 float2Int :: Float -> Int
 917 float2Int   (F# x) = I# (float2Int# x)
 918
 919 int2Float :: Int -> Float
 920 int2Float   (I# x) = F# (int2Float# x)
 921
 922 expFloat, logFloat, sqrtFloat :: Float -> Float
 923 sinFloat, cosFloat, tanFloat  :: Float -> Float
 924 asinFloat, acosFloat, atanFloat  :: Float -> Float
 925 sinhFloat, coshFloat, tanhFloat  :: Float -> Float
 926 expFloat    (F# x) = F# (expFloat# x)
 927 logFloat    (F# x) = F# (logFloat# x)
 928 sqrtFloat   (F# x) = F# (sqrtFloat# x)
 929 sinFloat    (F# x) = F# (sinFloat# x)
 930 cosFloat    (F# x) = F# (cosFloat# x)
 931 tanFloat    (F# x) = F# (tanFloat# x)
 932 asinFloat   (F# x) = F# (asinFloat# x)
 933 acosFloat   (F# x) = F# (acosFloat# x)
 934 atanFloat   (F# x) = F# (atanFloat# x)
 935 sinhFloat   (F# x) = F# (sinhFloat# x)
 936 coshFloat   (F# x) = F# (coshFloat# x)
 937 tanhFloat   (F# x) = F# (tanhFloat# x)
 938
 939 powerFloat :: Float -> Float -> Float
 940 powerFloat  (F# x) (F# y) = F# (powerFloat# x y)
 941
 942 -- definitions of the boxed PrimOps; these will be
 943 -- used in the case of partial applications, etc.
 944
 945 plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
 946 plusDouble   (D# x) (D# y) = D# (x +## y)
 947 minusDouble  (D# x) (D# y) = D# (x -## y)
 948 timesDouble  (D# x) (D# y) = D# (x *## y)
 949 divideDouble (D# x) (D# y) = D# (x /## y)
 950
 951 negateDouble :: Double -> Double
 952 negateDouble (D# x)        = D# (negateDouble# x)
 953
 954 gtDouble, geDouble, eqDouble, neDouble, leDouble, ltDouble :: Double -> Double -> Bool
 955 gtDouble    (D# x) (D# y) = x >## y
 956 geDouble    (D# x) (D# y) = x >=## y
 957 eqDouble    (D# x) (D# y) = x ==## y
 958 neDouble    (D# x) (D# y) = x /=## y
 959 ltDouble    (D# x) (D# y) = x <## y
 960 leDouble    (D# x) (D# y) = x <=## y
 961
 962 double2Int :: Double -> Int
 963 double2Int   (D# x) = I# (double2Int#   x)
 964
 965 int2Double :: Int -> Double
 966 int2Double   (I# x) = D# (int2Double#   x)
 967
 968 double2Float :: Double -> Float
 969 double2Float (D# x) = F# (double2Float# x)
 970 float2Double :: Float -> Double
 971 float2Double (F# x) = D# (float2Double# x)
 972
 973 expDouble, logDouble, sqrtDouble :: Double -> Double
 974 sinDouble, cosDouble, tanDouble  :: Double -> Double
 975 asinDouble, acosDouble, atanDouble  :: Double -> Double
 976 sinhDouble, coshDouble, tanhDouble  :: Double -> Double
 977 expDouble    (D# x) = D# (expDouble# x)
 978 logDouble    (D# x) = D# (logDouble# x)
 979 sqrtDouble   (D# x) = D# (sqrtDouble# x)
 980 sinDouble    (D# x) = D# (sinDouble# x)
 981 cosDouble    (D# x) = D# (cosDouble# x)
 982 tanDouble    (D# x) = D# (tanDouble# x)
 983 asinDouble   (D# x) = D# (asinDouble# x)
 984 acosDouble   (D# x) = D# (acosDouble# x)
 985 atanDouble   (D# x) = D# (atanDouble# x)
 986 sinhDouble   (D# x) = D# (sinhDouble# x)
 987 coshDouble   (D# x) = D# (coshDouble# x)
 988 tanhDouble   (D# x) = D# (tanhDouble# x)
 989
 990 powerDouble :: Double -> Double -> Double
 991 powerDouble  (D# x) (D# y) = D# (x **## y)
 992 \end{code}