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