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