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