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