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
  25 import PrelArr          ( Array, array, (!) )
  26 import PrelIOBase       ( unsafePerformIO )
  27 import Ix               ( Ix(..) )
  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 x    = 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# d#)
 297       = case decodeDouble# d#   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 x    = 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 fromRealFrac ::
 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 fromRealFrac    :: (RealFrac a, Fractional b) => a -> b
 358 fromRealFrac    =  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     toEnum         =  fromIntegral
 383     fromEnum       =  fromInteger . truncate   -- may overflow
 384     enumFrom       =  numericEnumFrom
 385     enumFromThen   =  numericEnumFromThen
 386     enumFromThenTo =  numericEnumFromThenTo
 387
 388 instance  Enum Double  where
 389     toEnum         =  fromIntegral
 390     fromEnum       =  fromInteger . truncate   -- may overflow
 391     enumFrom       =  numericEnumFrom
 392     enumFromThen   =  numericEnumFromThen
 393     enumFromThenTo =  numericEnumFromThenTo
 394
 395 numericEnumFrom         :: (Real a) => a -> [a]
 396 numericEnumFromThen     :: (Real a) => a -> a -> [a]
 397 numericEnumFromThenTo   :: (Real a) => a -> a -> a -> [a]
 398 numericEnumFrom         =  iterate (+1)
 399 numericEnumFromThen n m =  iterate (+(m-n)) n
 400 numericEnumFromThenTo n m p = takeWhile (if m >= n then (<= p) else (>= p))
 401                                       (numericEnumFromThen n m)
 402 \end{code}
 403
 404 @approxRational@, applied to two real fractional numbers x and epsilon,
 405 returns the simplest rational number within epsilon of x.  A rational
 406 number n%d in reduced form is said to be simpler than another n'%d' if
 407 abs n <= abs n' && d <= d'.  Any real interval contains a unique
 408 simplest rational; here, for simplicity, we assume a closed rational
 409 interval.  If such an interval includes at least one whole number, then
 410 the simplest rational is the absolutely least whole number.  Otherwise,
 411 the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
 412 and abs r' < d', and the simplest rational is q%1 + the reciprocal of
 413 the simplest rational between d'%r' and d%r.
 414
 415 \begin{code}
 416 approxRational          :: (RealFrac a) => a -> a -> Rational
 417 approxRational x eps    =  simplest (x-eps) (x+eps)
 418         where simplest x y | y < x      =  simplest y x
 419                            | x == y     =  xr
 420                            | x > 0      =  simplest' n d n' d'
 421                            | y < 0      =  - simplest' (-n') d' (-n) d
 422                            | otherwise  =  0 :% 1
 423                                         where xr  = toRational x
 424                                               n   = numerator xr
 425                                               d   = denominator xr
 426                                               nd' = toRational y
 427                                               n'  = numerator nd'
 428                                               d'  = denominator nd'
 429
 430               simplest' n d n' d'       -- assumes 0 < n%d < n'%d'
 431                         | r == 0     =  q :% 1
 432                         | q /= q'    =  (q+1) :% 1
 433                         | otherwise  =  (q*n''+d'') :% n''
 434                                      where (q,r)      =  quotRem n d
 435                                            (q',r')    =  quotRem n' d'
 436                                            nd''       =  simplest' d' r' d r
 437                                            n''        =  numerator nd''
 438                                            d''        =  denominator nd''
 439 \end{code}
 440
 441
 442 \begin{code}
 443 instance  (Integral a)  => Ord (Ratio a)  where
 444     (x:%y) <= (x':%y')  =  x * y' <= x' * y
 445     (x:%y) <  (x':%y')  =  x * y' <  x' * y
 446
 447 instance  (Integral a)  => Num (Ratio a)  where
 448     (x:%y) + (x':%y')   =  reduce (x*y' + x'*y) (y*y')
 449     (x:%y) - (x':%y')   =  reduce (x*y' - x'*y) (y*y')
 450     (x:%y) * (x':%y')   =  reduce (x * x') (y * y')
 451     negate (x:%y)       =  (-x) :% y
 452     abs (x:%y)          =  abs x :% y
 453     signum (x:%y)       =  signum x :% 1
 454     fromInteger x       =  fromInteger x :% 1
 455
 456 instance  (Integral a)  => Real (Ratio a)  where
 457     toRational (x:%y)   =  toInteger x :% toInteger y
 458
 459 instance  (Integral a)  => Fractional (Ratio a)  where
 460     (x:%y) / (x':%y')   =  (x*y') % (y*x')
 461     recip (x:%y)        =  if x < 0 then (-y) :% (-x) else y :% x
 462     fromRational (x:%y) =  fromInteger x :% fromInteger y
 463
 464 instance  (Integral a)  => RealFrac (Ratio a)  where
 465     properFraction (x:%y) = (fromIntegral q, r:%y)
 466                             where (q,r) = quotRem x y
 467
 468 instance  (Integral a)  => Enum (Ratio a)  where
 469     enumFrom            =  iterate ((+)1)
 470     enumFromThen n m    =  iterate ((+)(m-n)) n
 471     toEnum n            =  fromIntegral n :% 1
 472     fromEnum            =  fromInteger . truncate
 473
 474 ratio_prec :: Int
 475 ratio_prec = 7
 476
 477 instance  (Integral a)  => Show (Ratio a)  where
 478     showsPrec p (x:%y)  =  showParen (p > ratio_prec)
 479                                (shows x . showString " % " . shows y)
 480 \end{code}
 481
 482 @showRational@ converts a Rational to a string that looks like a
 483 floating point number, but without converting to any floating type
 484 (because of the possible overflow).
 485
 486 From/by Lennart, 94/09/26
 487
 488 \begin{code}
 489 showRational :: Int -> Rational -> String
 490 showRational n r =
 491     if r == 0 then
 492         "0.0"
 493     else
 494         let (r', e) = normalize r
 495         in  prR n r' e
 496
 497 startExpExp = 4 :: Int
 498
 499 -- make sure 1 <= r < 10
 500 normalize :: Rational -> (Rational, Int)
 501 normalize r = if r < 1 then
 502                   case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
 503               else
 504                   norm startExpExp r 0
 505         where norm :: Int -> Rational -> Int -> (Rational, Int)
 506               -- Invariant: r*10^e == original r
 507               norm 0  r e = (r, e)
 508               norm ee r e =
 509                 let n = 10^ee
 510                     tn = 10^n
 511                 in  if r >= tn then norm ee (r/tn) (e+n) else norm (ee-1) r e
 512
 513 drop0 "" = ""
 514 drop0 (c:cs) = c : reverse (dropWhile (=='0') (reverse cs))
 515
 516 prR :: Int -> Rational -> Int -> String
 517 prR n r e | r <  1  = prR n (r*10) (e-1)                -- final adjustment
 518 prR n r e | r >= 10 = prR n (r/10) (e+1)
 519 prR n r e0 =
 520         let s = show ((round (r * 10^n))::Integer)
 521             e = e0+1
 522         in  if e > 0 && e < 8 then
 523                 take e s ++ "." ++ drop0 (drop e s)
 524             else if e <= 0 && e > -3 then
 525                 "0." ++ take (-e) (repeat '0') ++ drop0 s
 526             else
 527                 head s : "."++ drop0 (tail s) ++ "e" ++ show e0
 528 \end{code}
 529
 530 [In response to a request for documentation of how fromRational works,
 531 Joe Fasel writes:] A quite reasonable request!  This code was added to
 532 the Prelude just before the 1.2 release, when Lennart, working with an
 533 early version of hbi, noticed that (read . show) was not the identity
 534 for floating-point numbers.  (There was a one-bit error about half the
 535 time.)  The original version of the conversion function was in fact
 536 simply a floating-point divide, as you suggest above. The new version
 537 is, I grant you, somewhat denser.
 538
 539 Unfortunately, Joe's code doesn't work!  Here's an example:
 540
 541 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
 542
 543 This program prints
 544         0.0000000000000000
 545 instead of
 546         1.8217369128763981e-300
 547
 548 Lennart's code follows, and it works...
 549
 550 \begin{pseudocode}
 551 {-# SPECIALISE fromRat ::
 552         Rational -> Double,
 553         Rational -> Float #-}
 554 fromRat :: (RealFloat a) => Rational -> a
 555 fromRat x = x'
 556         where x' = f e
 557
 558 --              If the exponent of the nearest floating-point number to x
 559 --              is e, then the significand is the integer nearest xb^(-e),
 560 --              where b is the floating-point radix.  We start with a good
 561 --              guess for e, and if it is correct, the exponent of the
 562 --              floating-point number we construct will again be e.  If
 563 --              not, one more iteration is needed.
 564
 565               f e   = if e' == e then y else f e'
 566                       where y      = encodeFloat (round (x * (1 % b)^^e)) e
 567                             (_,e') = decodeFloat y
 568               b     = floatRadix x'
 569
 570 --              We obtain a trial exponent by doing a floating-point
 571 --              division of x's numerator by its denominator.  The
 572 --              result of this division may not itself be the ultimate
 573 --              result, because of an accumulation of three rounding
 574 --              errors.
 575
 576               (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
 577                                         / fromInteger (denominator x))
 578 \end{pseudocode}
 579
 580 Now, here's Lennart's code.
 581
 582 \begin{code}
 583 --fromRat :: (RealFloat a) => Rational -> a
 584 fromRat x =
 585     if x == 0 then encodeFloat 0 0              -- Handle exceptional cases
 586     else if x < 0 then - fromRat' (-x)          -- first.
 587     else fromRat' x
 588
 589 -- Conversion process:
 590 -- Scale the rational number by the RealFloat base until
 591 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
 592 -- Then round the rational to an Integer and encode it with the exponent
 593 -- that we got from the scaling.
 594 -- To speed up the scaling process we compute the log2 of the number to get
 595 -- a first guess of the exponent.
 596
 597 fromRat' :: (RealFloat a) => Rational -> a
 598 fromRat' x = r
 599   where b = floatRadix r
 600         p = floatDigits r
 601         (minExp0, _) = floatRange r
 602         minExp = minExp0 - p            -- the real minimum exponent
 603         xMin = toRational (expt b (p-1))
 604         xMax = toRational (expt b p)
 605         p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
 606         f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
 607         (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
 608         r = encodeFloat (round x') p'
 609
 610 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
 611 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
 612 scaleRat b minExp xMin xMax p x =
 613     if p <= minExp then
 614         (x, p)
 615     else if x >= xMax then
 616         scaleRat b minExp xMin xMax (p+1) (x/b)
 617     else if x < xMin  then
 618         scaleRat b minExp xMin xMax (p-1) (x*b)
 619     else
 620         (x, p)
 621
 622 -- Exponentiation with a cache for the most common numbers.
 623 minExpt = 0::Int
 624 maxExpt = 1100::Int
 625 expt :: Integer -> Int -> Integer
 626 expt base n =
 627     if base == 2 && n >= minExpt && n <= maxExpt then
 628         expts!n
 629     else
 630         base^n
 631 expts :: Array Int Integer
 632 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
 633
 634 -- Compute the (floor of the) log of i in base b.
 635 -- Simplest way would be just divide i by b until it's smaller then b, but that would
 636 -- be very slow!  We are just slightly more clever.
 637 integerLogBase :: Integer -> Integer -> Int
 638 integerLogBase b i =
 639      if i < b then
 640         0
 641      else
 642         -- Try squaring the base first to cut down the number of divisions.
 643         let l = 2 * integerLogBase (b*b) i
 644             doDiv :: Integer -> Int -> Int
 645             doDiv i l = if i < b then l else doDiv (i `div` b) (l+1)
 646         in  doDiv (i `div` (b^l)) l
 647 \end{code}
 648
 649 %*********************************************************
 650 %*                                                      *
 651 \subsection{Printing out numbers}
 652 %*                                                      *
 653 %*********************************************************
 654
 655 \begin{code}
 656 --Exported from std library Numeric, defined here to
 657 --avoid mut. rec. between PrelNum and Numeric.
 658 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
 659 showSigned showPos p x = if x < 0 then showParen (p > 6)
 660                                                  (showChar '-' . showPos (-x))
 661                                   else showPos x
 662
 663 showFloat x  =  showString (formatRealFloat FFGeneric Nothing x)
 664
 665 -- These are the format types.  This type is not exported.
 666
 667 data FFFormat = FFExponent | FFFixed | FFGeneric --no need: deriving (Eq, Ord, Show)
 668
 669 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
 670 formatRealFloat fmt decs x = s
 671  where
 672   base = 10
 673   s = if isNaN x
 674       then "NaN"
 675       else
 676        if isInfinite x then
 677           if x < 0 then "-Infinity" else "Infinity"
 678        else
 679           if x < 0 || isNegativeZero x then
 680             '-':doFmt fmt (floatToDigits (toInteger base) (-x))
 681           else
 682             doFmt fmt (floatToDigits (toInteger base) x)
 683
 684   doFmt fmt (is, e) =
 685     let ds = map intToDigit is in
 686     case fmt of
 687      FFGeneric ->
 688       doFmt (if e <0 || e > 7 then FFExponent else FFFixed)
 689             (is,e)
 690      FFExponent ->
 691       case decs of
 692        Nothing ->
 693         let e' = if e==0 then 0 else e-1 in
 694         (case ds of
 695           [d]    -> d : ".0e"
 696           (d:ds) -> d : '.' : ds ++ "e") ++ show e'
 697        Just dec ->
 698         let dec' = max dec 1 in
 699         case is of
 700          [0] -> '0':'.':take dec' (repeat '0') ++ "e0"
 701          _ ->
 702           let
 703            (ei,is') = roundTo base (dec'+1) is
 704            d:ds = map intToDigit (if ei > 0 then init is' else is')
 705           in
 706           d:'.':ds ++ 'e':show (e-1+ei)
 707      FFFixed ->
 708       let
 709        mk0 ls = case ls of { "" -> "0" ; _ -> ls}
 710       in
 711       case decs of
 712        Nothing ->
 713          let
 714           f 0 s ds = mk0 (reverse s) ++ '.':mk0 ds
 715           f n s "" = f (n-1) ('0':s) ""
 716           f n s (d:ds) = f (n-1) (d:s) ds
 717          in
 718          f e "" ds
 719        Just dec ->
 720         let dec' = max dec 1 in
 721         if e >= 0 then
 722          let
 723           (ei,is') = roundTo base (dec' + e) is
 724           (ls,rs)  = splitAt (e+ei) (map intToDigit is')
 725          in
 726          mk0 ls ++ (if null rs then "" else '.':rs)
 727         else
 728          let
 729           (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
 730           d:ds = map intToDigit (if ei > 0 then is' else 0:is')
 731          in
 732          d : '.' : ds
 733
 734
 735 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
 736 roundTo base d is =
 737  let
 738   v = f d is
 739  in
 740  case v of
 741   (0,is) -> v
 742   (1,is) -> (1, 1:is)
 743  where
 744   b2 = base `div` 2
 745
 746   f n [] = (0, replicate n 0)
 747   f 0 (i:_) = (if i>=b2 then 1 else 0, [])
 748   f d (i:is) =
 749     let
 750      (c,ds) = f (d-1) is
 751      i' = c + i
 752     in
 753     if i' == base then (1,0:ds) else (0,i':ds)
 754
 755 --
 756 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
 757 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
 758 -- This version uses a much slower logarithm estimator. It should be improved.
 759
 760 -- This function returns a list of digits (Ints in [0..base-1]) and an
 761 -- exponent.
 762 --floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
 763 floatToDigits _ 0 = ([0], 0)
 764 floatToDigits base x =
 765  let
 766   (f0, e0) = decodeFloat x
 767   (minExp0, _) = floatRange x
 768   p = floatDigits x
 769   b = floatRadix x
 770   minExp = minExp0 - p -- the real minimum exponent
 771   -- Haskell requires that f be adjusted so denormalized numbers
 772   -- will have an impossibly low exponent.  Adjust for this.
 773   (f, e) =
 774    let n = minExp - e0 in
 775    if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
 776   (r, s, mUp, mDn) =
 777    if e >= 0 then
 778     let be = b^ e in
 779     if f == b^(p-1) then
 780       (f*be*b*2, 2*b, be*b, b)
 781     else
 782       (f*be*2, 2, be, be)
 783    else
 784     if e > minExp && f == b^(p-1) then
 785       (f*b*2, b^(-e+1)*2, b, 1)
 786     else
 787       (f*2, b^(-e)*2, 1, 1)
 788   k =
 789    let
 790     k0 =
 791      if b == 2 && base == 10 then
 792         -- logBase 10 2 is slightly bigger than 3/10 so
 793         -- the following will err on the low side.  Ignoring
 794         -- the fraction will make it err even more.
 795         -- Haskell promises that p-1 <= logBase b f < p.
 796         (p - 1 + e0) * 3 `div` 10
 797      else
 798         ceiling ((log (fromInteger (f+1)) +
 799                  fromInt e * log (fromInteger b)) /
 800                   fromInt e * log (fromInteger b))
 801
 802     fixup n =
 803       if n >= 0 then
 804         if r + mUp <= expt base n * s then n else fixup (n+1)
 805       else
 806         if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
 807    in
 808    fixup k0
 809
 810   gen ds rn sN mUpN mDnN =
 811    let
 812     (dn, rn') = (rn * base) `divMod` sN
 813     mUpN' = mUpN * base
 814     mDnN' = mDnN * base
 815    in
 816    case (rn' < mDnN', rn' + mUpN' > sN) of
 817     (True,  False) -> dn : ds
 818     (False, True)  -> dn+1 : ds
 819     (True,  True)  -> if rn' * 2 < sN then dn : ds else dn+1 : ds
 820     (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
 821
 822   rds =
 823    if k >= 0 then
 824       gen [] r (s * expt base k) mUp mDn
 825    else
 826      let bk = expt base (-k) in
 827      gen [] (r * bk) s (mUp * bk) (mDn * bk)
 828  in
 829  (map toInt (reverse rds), k)
 830
 831 \end{code}
 832
 833 %*********************************************************
 834 %*                                                      *
 835 \subsection{Numeric primops}
 836 %*                                                      *
 837 %*********************************************************
 838
 839 Definitions of the boxed PrimOps; these will be
 840 used in the case of partial applications, etc.
 841
 842 \begin{code}
 843 plusFloat   (F# x) (F# y) = F# (plusFloat# x y)
 844 minusFloat  (F# x) (F# y) = F# (minusFloat# x y)
 845 timesFloat  (F# x) (F# y) = F# (timesFloat# x y)
 846 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
 847 negateFloat (F# x)        = F# (negateFloat# x)
 848
 849 gtFloat     (F# x) (F# y) = gtFloat# x y
 850 geFloat     (F# x) (F# y) = geFloat# x y
 851 eqFloat     (F# x) (F# y) = eqFloat# x y
 852 neFloat     (F# x) (F# y) = neFloat# x y
 853 ltFloat     (F# x) (F# y) = ltFloat# x y
 854 leFloat     (F# x) (F# y) = leFloat# x y
 855
 856 float2Int   (F# x) = I# (float2Int# x)
 857 int2Float   (I# x) = F# (int2Float# x)
 858
 859 expFloat    (F# x) = F# (expFloat# x)
 860 logFloat    (F# x) = F# (logFloat# x)
 861 sqrtFloat   (F# x) = F# (sqrtFloat# x)
 862 sinFloat    (F# x) = F# (sinFloat# x)
 863 cosFloat    (F# x) = F# (cosFloat# x)
 864 tanFloat    (F# x) = F# (tanFloat# x)
 865 asinFloat   (F# x) = F# (asinFloat# x)
 866 acosFloat   (F# x) = F# (acosFloat# x)
 867 atanFloat   (F# x) = F# (atanFloat# x)
 868 sinhFloat   (F# x) = F# (sinhFloat# x)
 869 coshFloat   (F# x) = F# (coshFloat# x)
 870 tanhFloat   (F# x) = F# (tanhFloat# x)
 871
 872 powerFloat  (F# x) (F# y) = F# (powerFloat# x y)
 873
 874 -- definitions of the boxed PrimOps; these will be
 875 -- used in the case of partial applications, etc.
 876
 877 plusDouble   (D# x) (D# y) = D# (x +## y)
 878 minusDouble  (D# x) (D# y) = D# (x -## y)
 879 timesDouble  (D# x) (D# y) = D# (x *## y)
 880 divideDouble (D# x) (D# y) = D# (x /## y)
 881 negateDouble (D# x)        = D# (negateDouble# x)
 882
 883 gtDouble    (D# x) (D# y) = x >## y
 884 geDouble    (D# x) (D# y) = x >=## y
 885 eqDouble    (D# x) (D# y) = x ==## y
 886 neDouble    (D# x) (D# y) = x /=## y
 887 ltDouble    (D# x) (D# y) = x <## y
 888 leDouble    (D# x) (D# y) = x <=## y
 889
 890 double2Int   (D# x) = I# (double2Int#   x)
 891 int2Double   (I# x) = D# (int2Double#   x)
 892 double2Float (D# x) = F# (double2Float# x)
 893 float2Double (F# x) = D# (float2Double# x)
 894
 895 expDouble    (D# x) = D# (expDouble# x)
 896 logDouble    (D# x) = D# (logDouble# x)
 897 sqrtDouble   (D# x) = D# (sqrtDouble# x)
 898 sinDouble    (D# x) = D# (sinDouble# x)
 899 cosDouble    (D# x) = D# (cosDouble# x)
 900 tanDouble    (D# x) = D# (tanDouble# x)
 901 asinDouble   (D# x) = D# (asinDouble# x)
 902 acosDouble   (D# x) = D# (acosDouble# x)
 903 atanDouble   (D# x) = D# (atanDouble# x)
 904 sinhDouble   (D# x) = D# (sinhDouble# x)
 905 coshDouble   (D# x) = D# (coshDouble# x)
 906 tanhDouble   (D# x) = D# (tanhDouble# x)
 907
 908 powerDouble  (D# x) (D# y) = D# (x **## y)
 909 \end{code}