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