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