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