ghc/lib/std/PrelNum.lhs

   1 %
   2 % (c) The AQUA Project, Glasgow University, 1994-1996
   3 %
   4
   5 \section[PrelNum]{Module @PrelNum@}
   6
   7 Numeric part of the prelude.
   8
   9 It's rather big!
  10
  11 \begin{code}
  12 {-# OPTIONS -fno-implicit-prelude -#include "cbits/floatExtreme.h" #-}
  13 {-# OPTIONS -H20m #-}
  14
  15 #include "../includes/ieee-flpt.h"
  16
  17 \end{code}
  18
  19 \begin{code}
  20 module PrelNum where
  21
  22 import PrelBase
  23 import PrelGHC
  24 import {-# SOURCE #-} PrelErr ( error )
  25 import PrelList
  26 import PrelMaybe
  27
  28 import PrelArr          ( Array, array, (!) )
  29 import PrelIOBase       ( unsafePerformIO )
  30 import Ix               ( Ix(..) )
  31 import PrelCCall        ()      -- we need the definitions of CCallable and
  32                                 -- CReturnable for the _ccall_s herein.
  33
  34
  35 infixr 8  ^, ^^, **
  36 infixl 7  /, %, `quot`, `rem`, `div`, `mod`
  37 \end{code}
  38
  39
  40 %*********************************************************
  41 %*                                                      *
  42 \subsection{Standard numeric classes}
  43 %*                                                      *
  44 %*********************************************************
  45
  46 \begin{code}
  47 class  (Num a, Ord a) => Real a  where
  48     toRational          ::  a -> Rational
  49
  50 class  (Real a, Enum a) => Integral a  where
  51     quot, rem, div, mod :: a -> a -> a
  52     quotRem, divMod     :: a -> a -> (a,a)
  53     toInteger           :: a -> Integer
  54     toInt               :: a -> Int -- partain: Glasgow extension
  55
  56     n `quot` d          =  q  where (q,r) = quotRem n d
  57     n `rem` d           =  r  where (q,r) = quotRem n d
  58     n `div` d           =  q  where (q,r) = divMod n d
  59     n `mod` d           =  r  where (q,r) = divMod n d
  60     divMod n d          =  if signum r == negate (signum d) then (q-1, r+d) else qr
  61                            where qr@(q,r) = quotRem n d
  62
  63 class  (Num a) => Fractional a  where
  64     (/)                 :: a -> a -> a
  65     recip               :: a -> a
  66     fromRational        :: Rational -> a
  67
  68     recip x             =  1 / x
  69
  70 class  (Fractional a) => Floating a  where
  71     pi                  :: a
  72     exp, log, sqrt      :: a -> a
  73     (**), logBase       :: a -> a -> a
  74     sin, cos, tan       :: a -> a
  75     asin, acos, atan    :: a -> a
  76     sinh, cosh, tanh    :: a -> a
  77     asinh, acosh, atanh :: a -> a
  78
  79     x ** y              =  exp (log x * y)
  80     logBase x y         =  log y / log x
  81     sqrt x              =  x ** 0.5
  82     tan  x              =  sin  x / cos  x
  83     tanh x              =  sinh x / cosh x
  84
  85 class  (Real a, Fractional a) => RealFrac a  where
  86     properFraction      :: (Integral b) => a -> (b,a)
  87     truncate, round     :: (Integral b) => a -> b
  88     ceiling, floor      :: (Integral b) => a -> b
  89
  90     truncate x          =  m  where (m,_) = properFraction x
  91
  92     round x             =  let (n,r) = properFraction x
  93                                m     = if r < 0 then n - 1 else n + 1
  94                            in case signum (abs r - 0.5) of
  95                                 -1 -> n
  96                                 0  -> if even n then n else m
  97                                 1  -> m
  98
  99     ceiling x           =  if r > 0 then n + 1 else n
 100                            where (n,r) = properFraction x
 101
 102     floor x             =  if r < 0 then n - 1 else n
 103                            where (n,r) = properFraction x
 104
 105 class  (RealFrac a, Floating a) => RealFloat a  where
 106     floatRadix          :: a -> Integer
 107     floatDigits         :: a -> Int
 108     floatRange          :: a -> (Int,Int)
 109     decodeFloat         :: a -> (Integer,Int)
 110     encodeFloat         :: Integer -> Int -> a
 111     exponent            :: a -> Int
 112     significand         :: a -> a
 113     scaleFloat          :: Int -> a -> a
 114     isNaN, isInfinite, isDenormalized, isNegativeZero, isIEEE
 115                         :: a -> Bool
 116
 117     exponent x          =  if m == 0 then 0 else n + floatDigits x
 118                            where (m,n) = decodeFloat x
 119
 120     significand x       =  encodeFloat m (negate (floatDigits x))
 121                            where (m,_) = decodeFloat x
 122
 123     scaleFloat k x      =  encodeFloat m (n+k)
 124                            where (m,n) = decodeFloat x
 125 \end{code}
 126
 127 %*********************************************************
 128 %*                                                      *
 129 \subsection{Overloaded numeric functions}
 130 %*                                                      *
 131 %*********************************************************
 132
 133 \begin{code}
 134 even, odd       :: (Integral a) => a -> Bool
 135 even n          =  n `rem` 2 == 0
 136 odd             =  not . even
 137
 138 {-# SPECIALISE gcd ::
 139         Int -> Int -> Int,
 140         Integer -> Integer -> Integer #-}
 141 gcd             :: (Integral a) => a -> a -> a
 142 gcd 0 0         =  error "Prelude.gcd: gcd 0 0 is undefined"
 143 gcd x y         =  gcd' (abs x) (abs y)
 144                    where gcd' x 0  =  x
 145                          gcd' x y  =  gcd' y (x `rem` y)
 146
 147 {-# SPECIALISE lcm ::
 148         Int -> Int -> Int,
 149         Integer -> Integer -> Integer #-}
 150 lcm             :: (Integral a) => a -> a -> a
 151 lcm _ 0         =  0
 152 lcm 0 _         =  0
 153 lcm x y         =  abs ((x `quot` (gcd x y)) * y)
 154
 155 {-# SPECIALISE (^) ::
 156         Integer -> Integer -> Integer,
 157         Integer -> Int -> Integer,
 158         Int -> Int -> Int #-}
 159 (^)             :: (Num a, Integral b) => a -> b -> a
 160 x ^ 0           =  1
 161 x ^ n | n > 0   =  f x (n-1) x
 162                    where f _ 0 y = y
 163                          f x n y = g x n  where
 164                                    g x n | even n  = g (x*x) (n `quot` 2)
 165                                          | otherwise = f x (n-1) (x*y)
 166 _ ^ _           = error "Prelude.^: negative exponent"
 167
 168 {-# SPECIALISE (^^) ::
 169         Double -> Int -> Double,
 170         Rational -> Int -> Rational #-}
 171 (^^)            :: (Fractional a, Integral b) => a -> b -> a
 172 x ^^ n          =  if n >= 0 then x^n else recip (x^(negate n))
 173
 174 {-# SPECIALIZE fromIntegral ::
 175     Int         -> Rational,
 176     Integer     -> Rational,
 177     Int         -> Int,
 178     Int         -> Integer,
 179     Int         -> Float,
 180     Int         -> Double,
 181     Integer     -> Int,
 182     Integer     -> Integer,
 183     Integer     -> Float,
 184     Integer     -> Double #-}
 185 fromIntegral    :: (Integral a, Num b) => a -> b
 186 fromIntegral    =  fromInteger . toInteger
 187
 188 {-# SPECIALIZE fromRealFrac ::
 189     Double      -> Rational,
 190     Rational    -> Double,
 191     Float       -> Rational,
 192     Rational    -> Float,
 193     Rational    -> Rational,
 194     Double      -> Double,
 195     Double      -> Float,
 196     Float       -> Float,
 197     Float       -> Double #-}
 198 fromRealFrac    :: (RealFrac a, Fractional b) => a -> b
 199 fromRealFrac    =  fromRational . toRational
 200
 201 atan2           :: (RealFloat a) => a -> a -> a
 202 atan2 y x       =  case (signum y, signum x) of
 203                         ( 0, 1) ->  0
 204                         ( 1, 0) ->  pi/2
 205                         ( 0,-1) ->  pi
 206                         (-1, 0) ->  (negate pi)/2
 207                         ( _, 1) ->  atan (y/x)
 208                         ( _,-1) ->  atan (y/x) + pi
 209                         ( 0, 0) ->  error "Prelude.atan2: atan2 of origin"
 210 \end{code}
 211
 212
 213 %*********************************************************
 214 %*                                                      *
 215 \subsection{Instances for @Int@}
 216 %*                                                      *
 217 %*********************************************************
 218
 219 \begin{code}
 220 instance  Real Int  where
 221     toRational x        =  toInteger x % 1
 222
 223 instance  Integral Int  where
 224     a@(I# _) `quotRem` b@(I# _) = (a `quotInt` b, a `remInt` b)
 225     -- OK, so I made it a little stricter.  Shoot me.  (WDP 94/10)
 226
 227     -- Following chks for zero divisor are non-standard (WDP)
 228     a `quot` b  =  if b /= 0
 229                    then a `quotInt` b
 230                    else error "Integral.Int.quot{PreludeCore}: divide by 0\n"
 231     a `rem` b   =  if b /= 0
 232                    then a `remInt` b
 233                    else error "Integral.Int.rem{PreludeCore}: divide by 0\n"
 234
 235     x `div` y = if x > 0 && y < 0       then quotInt (x-y-1) y
 236                 else if x < 0 && y > 0  then quotInt (x-y+1) y
 237                 else quotInt x y
 238     x `mod` y = if x > 0 && y < 0 || x < 0 && y > 0 then
 239                     if r/=0 then r+y else 0
 240                 else
 241                     r
 242               where r = remInt x y
 243
 244     divMod x@(I# _) y@(I# _) = (x `div` y, x `mod` y)
 245     -- Stricter.  Sorry if you don't like it.  (WDP 94/10)
 246
 247 --OLD:   even x = eqInt (x `mod` 2) 0
 248 --OLD:   odd x  = neInt (x `mod` 2) 0
 249
 250     toInteger (I# n#) = int2Integer# n#  -- give back a full-blown Integer
 251     toInt x           = x
 252
 253 \end{code}
 254
 255
 256 %*********************************************************
 257 %*                                                      *
 258 \subsection{Type @Integer@}
 259 %*                                                      *
 260 %*********************************************************
 261
 262 These types are used to return from integer primops
 263
 264 \begin{code}
 265 data Return2GMPs     = Return2GMPs     Int# Int# ByteArray# Int# Int# ByteArray#
 266 data ReturnIntAndGMP = ReturnIntAndGMP Int# Int# Int# ByteArray#
 267 \end{code}
 268
 269 Instances
 270
 271 \begin{code}
 272 instance  Eq Integer  where
 273     (J# a1 s1 d1) == (J# a2 s2 d2)
 274       = (cmpInteger# a1 s1 d1 a2 s2 d2) ==# 0#
 275
 276     (J# a1 s1 d1) /= (J# a2 s2 d2)
 277       = (cmpInteger# a1 s1 d1 a2 s2 d2) /=# 0#
 278
 279 instance  Ord Integer  where
 280     (J# a1 s1 d1) <= (J# a2 s2 d2)
 281       = (cmpInteger# a1 s1 d1 a2 s2 d2) <=# 0#
 282
 283     (J# a1 s1 d1) <  (J# a2 s2 d2)
 284       = (cmpInteger# a1 s1 d1 a2 s2 d2) <# 0#
 285
 286     (J# a1 s1 d1) >= (J# a2 s2 d2)
 287       = (cmpInteger# a1 s1 d1 a2 s2 d2) >=# 0#
 288
 289     (J# a1 s1 d1) >  (J# a2 s2 d2)
 290       = (cmpInteger# a1 s1 d1 a2 s2 d2) ># 0#
 291
 292     x@(J# a1 s1 d1) `max` y@(J# a2 s2 d2)
 293       = if ((cmpInteger# a1 s1 d1 a2 s2 d2) ># 0#) then x else y
 294
 295     x@(J# a1 s1 d1) `min` y@(J# a2 s2 d2)
 296       = if ((cmpInteger# a1 s1 d1 a2 s2 d2) <# 0#) then x else y
 297
 298     compare (J# a1 s1 d1) (J# a2 s2 d2)
 299        = case cmpInteger# a1 s1 d1 a2 s2 d2 of { res# ->
 300          if res# <# 0# then LT else
 301          if res# ># 0# then GT else EQ
 302          }
 303
 304 instance  Num Integer  where
 305     (+) (J# a1 s1 d1) (J# a2 s2 d2)
 306       = plusInteger# a1 s1 d1 a2 s2 d2
 307
 308     (-) (J# a1 s1 d1) (J# a2 s2 d2)
 309       = minusInteger# a1 s1 d1 a2 s2 d2
 310
 311     negate (J# a s d) = negateInteger# a s d
 312
 313     (*) (J# a1 s1 d1) (J# a2 s2 d2)
 314       = timesInteger# a1 s1 d1 a2 s2 d2
 315
 316     -- ORIG: abs n = if n >= 0 then n else -n
 317
 318     abs n@(J# a1 s1 d1)
 319       = case 0 of { J# a2 s2 d2 ->
 320         if (cmpInteger# a1 s1 d1 a2 s2 d2) >=# 0#
 321         then n
 322         else negateInteger# a1 s1 d1
 323         }
 324
 325     signum n@(J# a1 s1 d1)
 326       = case 0 of { J# a2 s2 d2 ->
 327         let
 328             cmp = cmpInteger# a1 s1 d1 a2 s2 d2
 329         in
 330         if      cmp >#  0# then 1
 331         else if cmp ==# 0# then 0
 332         else                    (negate 1)
 333         }
 334
 335     fromInteger x       =  x
 336
 337     fromInt (I# n#)     =  int2Integer# n# -- gives back a full-blown Integer
 338
 339 instance  Real Integer  where
 340     toRational x        =  x % 1
 341
 342 instance  Integral Integer where
 343     quotRem (J# a1 s1 d1) (J# a2 s2 d2)
 344       = case (quotRemInteger# a1 s1 d1 a2 s2 d2) of
 345           Return2GMPs a3 s3 d3 a4 s4 d4
 346             -> (J# a3 s3 d3, J# a4 s4 d4)
 347
 348 {- USING THE UNDERLYING "GMP" CODE IS DUBIOUS FOR NOW:
 349
 350     divMod (J# a1 s1 d1) (J# a2 s2 d2)
 351       = case (divModInteger# a1 s1 d1 a2 s2 d2) of
 352           Return2GMPs a3 s3 d3 a4 s4 d4
 353             -> (J# a3 s3 d3, J# a4 s4 d4)
 354 -}
 355     toInteger n      = n
 356     toInt (J# a s d) = case (integer2Int# a s d) of { n# -> I# n# }
 357
 358     -- the rest are identical to the report default methods;
 359     -- you get slightly better code if you let the compiler
 360     -- see them right here:
 361     n `quot` d  =  if d /= 0 then q else
 362                      error "Integral.Integer.quot{PreludeCore}: divide by 0\n"
 363                    where (q,r) = quotRem n d
 364     n `rem` d   =  if d /= 0 then r else
 365                      error "Integral.Integer.quot{PreludeCore}: divide by 0\n"
 366                    where (q,r) = quotRem n d
 367     n `div` d   =  q  where (q,r) = divMod n d
 368     n `mod` d   =  r  where (q,r) = divMod n d
 369
 370     divMod n d  =  case (quotRem n d) of { qr@(q,r) ->
 371                    if signum r == negate (signum d) then (q - 1, r+d) else qr }
 372                    -- Case-ified by WDP 94/10
 373
 374 instance  Enum Integer  where
 375     toEnum n             =  toInteger n
 376     fromEnum n           =  toInt n
 377     enumFrom n           =  n : enumFrom (n + 1)
 378     enumFromThen m n     =  en' m (n - m)
 379                             where en' m n = m : en' (m + n) n
 380     enumFromTo n m       =  takeWhile (<= m) (enumFrom n)
 381     enumFromThenTo n m p =  takeWhile (if m >= n then (<= p) else (>= p))
 382                                       (enumFromThen n m)
 383
 384 instance  Show Integer  where
 385     showsPrec   x = showSignedInteger x
 386     showList = showList__ (showsPrec 0)
 387
 388 instance  Ix Integer  where
 389     range (m,n)         =  [m..n]
 390     index b@(m,n) i
 391         | inRange b i   =  fromInteger (i - m)
 392         | otherwise     =  error "Integer.index: Index out of range."
 393     inRange (m,n) i     =  m <= i && i <= n
 394
 395 integer_0, integer_1, integer_2, integer_m1 :: Integer
 396 integer_0  = int2Integer# 0#
 397 integer_1  = int2Integer# 1#
 398 integer_2  = int2Integer# 2#
 399 integer_m1 = int2Integer# (negateInt# 1#)
 400 \end{code}
 401
 402
 403 %*********************************************************
 404 %*                                                      *
 405 \subsection{Type @Float@}
 406 %*                                                      *
 407 %*********************************************************
 408
 409 \begin{code}
 410 instance Eq Float where
 411     (F# x) == (F# y) = x `eqFloat#` y
 412
 413 instance Ord Float where
 414     (F# x) `compare` (F# y) | x `ltFloat#` y = LT
 415                             | x `eqFloat#` y = EQ
 416                             | otherwise      = GT
 417
 418     (F# x) <  (F# y) = x `ltFloat#`  y
 419     (F# x) <= (F# y) = x `leFloat#`  y
 420     (F# x) >= (F# y) = x `geFloat#`  y
 421     (F# x) >  (F# y) = x `gtFloat#`  y
 422
 423 instance  Num Float  where
 424     (+)         x y     =  plusFloat x y
 425     (-)         x y     =  minusFloat x y
 426     negate      x       =  negateFloat x
 427     (*)         x y     =  timesFloat x y
 428     abs x | x >= 0.0    =  x
 429           | otherwise   =  negateFloat x
 430     signum x | x == 0.0  = 0
 431              | x > 0.0   = 1
 432              | otherwise = negate 1
 433     fromInteger n       =  encodeFloat n 0
 434     fromInt i           =  int2Float i
 435
 436 instance  Real Float  where
 437     toRational x        =  (m%1)*(b%1)^^n
 438                            where (m,n) = decodeFloat x
 439                                  b     = floatRadix  x
 440
 441 instance  Fractional Float  where
 442     (/) x y             =  divideFloat x y
 443     fromRational x      =  fromRat x
 444     recip x             =  1.0 / x
 445
 446 instance  Floating Float  where
 447     pi                  =  3.141592653589793238
 448     exp x               =  expFloat x
 449     log x               =  logFloat x
 450     sqrt x              =  sqrtFloat x
 451     sin x               =  sinFloat x
 452     cos x               =  cosFloat x
 453     tan x               =  tanFloat x
 454     asin x              =  asinFloat x
 455     acos x              =  acosFloat x
 456     atan x              =  atanFloat x
 457     sinh x              =  sinhFloat x
 458     cosh x              =  coshFloat x
 459     tanh x              =  tanhFloat x
 460     (**) x y            =  powerFloat x y
 461     logBase x y         =  log y / log x
 462
 463     asinh x = log (x + sqrt (1.0+x*x))
 464     acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
 465     atanh x = log ((x+1.0) / sqrt (1.0-x*x))
 466
 467 instance  RealFrac Float  where
 468
 469     {-# SPECIALIZE properFraction :: Float -> (Int, Float) #-}
 470     {-# SPECIALIZE truncate :: Float -> Int #-}
 471     {-# SPECIALIZE round    :: Float -> Int #-}
 472     {-# SPECIALIZE ceiling  :: Float -> Int #-}
 473     {-# SPECIALIZE floor    :: Float -> Int #-}
 474
 475     {-# SPECIALIZE properFraction :: Float -> (Integer, Float) #-}
 476     {-# SPECIALIZE truncate :: Float -> Integer #-}
 477     {-# SPECIALIZE round    :: Float -> Integer #-}
 478     {-# SPECIALIZE ceiling  :: Float -> Integer #-}
 479     {-# SPECIALIZE floor    :: Float -> Integer #-}
 480
 481     properFraction x
 482       = case (decodeFloat x)      of { (m,n) ->
 483         let  b = floatRadix x     in
 484         if n >= 0 then
 485             (fromInteger m * fromInteger b ^ n, 0.0)
 486         else
 487             case (quotRem m (b^(negate n))) of { (w,r) ->
 488             (fromInteger w, encodeFloat r n)
 489             }
 490         }
 491
 492     truncate x  = case properFraction x of
 493                      (n,_) -> n
 494
 495     round x     = case properFraction x of
 496                      (n,r) -> let
 497                                 m         = if r < 0.0 then n - 1 else n + 1
 498                                 half_down = abs r - 0.5
 499                               in
 500                               case (compare half_down 0.0) of
 501                                 LT -> n
 502                                 EQ -> if even n then n else m
 503                                 GT -> m
 504
 505     ceiling x   = case properFraction x of
 506                     (n,r) -> if r > 0.0 then n + 1 else n
 507
 508     floor x     = case properFraction x of
 509                     (n,r) -> if r < 0.0 then n - 1 else n
 510
 511 instance  RealFloat Float  where
 512     floatRadix _        =  FLT_RADIX        -- from float.h
 513     floatDigits _       =  FLT_MANT_DIG     -- ditto
 514     floatRange _        =  (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
 515
 516     decodeFloat (F# f#)
 517       = case decodeFloat# f#    of
 518           ReturnIntAndGMP exp# a# s# d# ->
 519             (J# a# s# d#, I# exp#)
 520
 521     encodeFloat (J# a# s# d#) (I# e#)
 522       = case encodeFloat# a# s# d# e# of { flt# -> F# flt# }
 523
 524     exponent x          = case decodeFloat x of
 525                             (m,n) -> if m == 0 then 0 else n + floatDigits x
 526
 527     significand x       = case decodeFloat x of
 528                             (m,_) -> encodeFloat m (negate (floatDigits x))
 529
 530     scaleFloat k x      = case decodeFloat x of
 531                             (m,n) -> encodeFloat m (n+k)
 532     isNaN x =
 533       (0::Int) /= unsafePerformIO (_ccall_ isFloatNaN x) {- a _pure_function! -}
 534     isInfinite x =
 535       (0::Int) /= unsafePerformIO (_ccall_ isFloatInfinite x) {- ditto! -}
 536     isDenormalized x =
 537       (0::Int) /= unsafePerformIO (_ccall_ isFloatDenormalized x) -- ..
 538     isNegativeZero x =
 539       (0::Int) /= unsafePerformIO (_ccall_ isFloatNegativeZero x) -- ...
 540     isIEEE x    = True
 541
 542 instance  Show Float  where
 543     showsPrec   x = showSigned showFloat x
 544     showList = showList__ (showsPrec 0)
 545 \end{code}
 546
 547 %*********************************************************
 548 %*                                                      *
 549 \subsection{Type @Double@}
 550 %*                                                      *
 551 %*********************************************************
 552
 553 \begin{code}
 554 instance Eq Double where
 555     (D# x) == (D# y) = x ==## y
 556
 557 instance Ord Double where
 558     (D# x) `compare` (D# y) | x <## y   = LT
 559                             | x ==## y  = EQ
 560                             | otherwise = GT
 561
 562     (D# x) <  (D# y) = x <##  y
 563     (D# x) <= (D# y) = x <=## y
 564     (D# x) >= (D# y) = x >=## y
 565     (D# x) >  (D# y) = x >##  y
 566
 567 instance  Num Double  where
 568     (+)         x y     =  plusDouble x y
 569     (-)         x y     =  minusDouble x y
 570     negate      x       =  negateDouble x
 571     (*)         x y     =  timesDouble x y
 572     abs x | x >= 0.0    =  x
 573           | otherwise   =  negateDouble x
 574     signum x | x == 0.0  = 0
 575              | x > 0.0   = 1
 576              | otherwise = negate 1
 577     fromInteger n       =  encodeFloat n 0
 578     fromInt (I# n#)     =  case (int2Double# n#) of { d# -> D# d# }
 579
 580 instance  Real Double  where
 581     toRational x        =  (m%1)*(b%1)^^n
 582                            where (m,n) = decodeFloat x
 583                                  b     = floatRadix  x
 584
 585 instance  Fractional Double  where
 586     (/) x y             =  divideDouble x y
 587     fromRational x      =  fromRat x
 588     recip x             =  1.0 / x
 589
 590 instance  Floating Double  where
 591     pi                  =  3.141592653589793238
 592     exp x               =  expDouble x
 593     log x               =  logDouble x
 594     sqrt x              =  sqrtDouble x
 595     sin  x              =  sinDouble x
 596     cos  x              =  cosDouble x
 597     tan  x              =  tanDouble x
 598     asin x              =  asinDouble x
 599     acos x              =  acosDouble x
 600     atan x              =  atanDouble x
 601     sinh x              =  sinhDouble x
 602     cosh x              =  coshDouble x
 603     tanh x              =  tanhDouble x
 604     (**) x y            =  powerDouble x y
 605     logBase x y         =  log y / log x
 606
 607     asinh x = log (x + sqrt (1.0+x*x))
 608     acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
 609     atanh x = log ((x+1.0) / sqrt (1.0-x*x))
 610
 611 instance  RealFrac Double  where
 612
 613     {-# SPECIALIZE properFraction :: Double -> (Int, Double) #-}
 614     {-# SPECIALIZE truncate :: Double -> Int #-}
 615     {-# SPECIALIZE round    :: Double -> Int #-}
 616     {-# SPECIALIZE ceiling  :: Double -> Int #-}
 617     {-# SPECIALIZE floor    :: Double -> Int #-}
 618
 619     {-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
 620     {-# SPECIALIZE truncate :: Double -> Integer #-}
 621     {-# SPECIALIZE round    :: Double -> Integer #-}
 622     {-# SPECIALIZE ceiling  :: Double -> Integer #-}
 623     {-# SPECIALIZE floor    :: Double -> Integer #-}
 624
 625 #if defined(__UNBOXED_INSTANCES__)
 626     {-# SPECIALIZE properFraction :: Double -> (Int#, Double) #-}
 627     {-# SPECIALIZE truncate :: Double -> Int# #-}
 628     {-# SPECIALIZE round    :: Double -> Int# #-}
 629     {-# SPECIALIZE ceiling  :: Double -> Int# #-}
 630     {-# SPECIALIZE floor    :: Double -> Int# #-}
 631 #endif
 632
 633     properFraction x
 634       = case (decodeFloat x)      of { (m,n) ->
 635         let  b = floatRadix x     in
 636         if n >= 0 then
 637             (fromInteger m * fromInteger b ^ n, 0.0)
 638         else
 639             case (quotRem m (b^(negate n))) of { (w,r) ->
 640             (fromInteger w, encodeFloat r n)
 641             }
 642         }
 643
 644     truncate x  = case properFraction x of
 645                      (n,_) -> n
 646
 647     round x     = case properFraction x of
 648                      (n,r) -> let
 649                                 m         = if r < 0.0 then n - 1 else n + 1
 650                                 half_down = abs r - 0.5
 651                               in
 652                               case (compare half_down 0.0) of
 653                                 LT -> n
 654                                 EQ -> if even n then n else m
 655                                 GT -> m
 656
 657     ceiling x   = case properFraction x of
 658                     (n,r) -> if r > 0.0 then n + 1 else n
 659
 660     floor x     = case properFraction x of
 661                     (n,r) -> if r < 0.0 then n - 1 else n
 662
 663 instance  RealFloat Double  where
 664     floatRadix _        =  FLT_RADIX        -- from float.h
 665     floatDigits _       =  DBL_MANT_DIG     -- ditto
 666     floatRange _        =  (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
 667
 668     decodeFloat (D# d#)
 669       = case decodeDouble# d#   of
 670           ReturnIntAndGMP exp# a# s# d# ->
 671             (J# a# s# d#, I# exp#)
 672
 673     encodeFloat (J# a# s# d#) (I# e#)
 674       = case encodeDouble# a# s# d# e#  of { dbl# -> D# dbl# }
 675
 676     exponent x          = case decodeFloat x of
 677                             (m,n) -> if m == 0 then 0 else n + floatDigits x
 678
 679     significand x       = case decodeFloat x of
 680                             (m,_) -> encodeFloat m (negate (floatDigits x))
 681
 682     scaleFloat k x      = case decodeFloat x of
 683                             (m,n) -> encodeFloat m (n+k)
 684     isNaN x =
 685       (0::Int) /= unsafePerformIO (_ccall_ isDoubleNaN x) {- a _pure_function! -}
 686     isInfinite x =
 687       (0::Int) /= unsafePerformIO (_ccall_ isDoubleInfinite x) {- ditto -}
 688     isDenormalized x =
 689       (0::Int) /= unsafePerformIO (_ccall_ isDoubleDenormalized x) -- ..
 690     isNegativeZero x =
 691       (0::Int) /= unsafePerformIO (_ccall_ isDoubleNegativeZero x) -- ...
 692     isIEEE x    = True
 693
 694 instance  Show Double  where
 695     showsPrec   x = showSigned showFloat x
 696     showList = showList__ (showsPrec 0)
 697 \end{code}
 698
 699
 700 %*********************************************************
 701 %*                                                      *
 702 \subsection{Common code for @Float@ and @Double@}
 703 %*                                                      *
 704 %*********************************************************
 705
 706 The @Enum@ instances for Floats and Doubles are slightly unusual.
 707 The @toEnum@ function truncates numbers to Int.  The definitions
 708 of @enumFrom@ and @enumFromThen@ allow floats to be used in arithmetic
 709 series: [0,0.1 .. 1.0].  However, roundoff errors make these somewhat
 710 dubious.  This example may have either 10 or 11 elements, depending on
 711 how 0.1 is represented.
 712
 713 NOTE: The instances for Float and Double do not make use of the default
 714 methods for @enumFromTo@ and @enumFromThenTo@, as these rely on there being
 715 a `non-lossy' conversion to and from Ints. Instead we make use of the
 716 1.2 default methods (back in the days when Enum had Ord as a superclass)
 717 for these (@numericEnumFromTo@ and @numericEnumFromThenTo@ below.)
 718
 719 \begin{code}
 720 instance  Enum Float  where
 721     toEnum         =  fromIntegral
 722     fromEnum       =  fromInteger . truncate   -- may overflow
 723     enumFrom       =  numericEnumFrom
 724     enumFromThen   =  numericEnumFromThen
 725     enumFromThenTo =  numericEnumFromThenTo
 726
 727 instance  Enum Double  where
 728     toEnum         =  fromIntegral
 729     fromEnum       =  fromInteger . truncate   -- may overflow
 730     enumFrom       =  numericEnumFrom
 731     enumFromThen   =  numericEnumFromThen
 732     enumFromThenTo =  numericEnumFromThenTo
 733
 734 numericEnumFrom         :: (Real a) => a -> [a]
 735 numericEnumFromThen     :: (Real a) => a -> a -> [a]
 736 numericEnumFromThenTo   :: (Real a) => a -> a -> a -> [a]
 737 numericEnumFrom         =  iterate (+1)
 738 numericEnumFromThen n m =  iterate (+(m-n)) n
 739 numericEnumFromThenTo n m p = takeWhile (if m >= n then (<= p) else (>= p))
 740                                       (numericEnumFromThen n m)
 741 \end{code}
 742
 743
 744 %*********************************************************
 745 %*                                                      *
 746 \subsection{The @Ratio@ and @Rational@ types}
 747 %*                                                      *
 748 %*********************************************************
 749
 750 \begin{code}
 751 data  (Eval a, Integral a)      => Ratio a = !a :% !a  deriving (Eq)
 752 type  Rational          =  Ratio Integer
 753 \end{code}
 754
 755 \begin{code}
 756 {-# SPECIALISE (%) :: Integer -> Integer -> Rational #-}
 757
 758 (%)                     :: (Integral a) => a -> a -> Ratio a
 759 numerator, denominator  :: (Integral a) => Ratio a -> a
 760 approxRational          :: (RealFrac a) => a -> a -> Rational
 761
 762 \end{code}
 763
 764 \tr{reduce} is a subsidiary function used only in this module .
 765 It normalises a ratio by dividing both numerator and denominator by
 766 their greatest common divisor.
 767
 768 \begin{code}
 769 reduce x 0              =  error "{Ratio.%}: zero denominator"
 770 reduce x y              =  (x `quot` d) :% (y `quot` d)
 771                            where d = gcd x y
 772 \end{code}
 773
 774 \begin{code}
 775 x % y                   =  reduce (x * signum y) (abs y)
 776
 777 numerator (x:%y)        =  x
 778
 779 denominator (x:%y)      =  y
 780 \end{code}
 781
 782
 783 @approxRational@, applied to two real fractional numbers x and epsilon,
 784 returns the simplest rational number within epsilon of x.  A rational
 785 number n%d in reduced form is said to be simpler than another n'%d' if
 786 abs n <= abs n' && d <= d'.  Any real interval contains a unique
 787 simplest rational; here, for simplicity, we assume a closed rational
 788 interval.  If such an interval includes at least one whole number, then
 789 the simplest rational is the absolutely least whole number.  Otherwise,
 790 the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
 791 and abs r' < d', and the simplest rational is q%1 + the reciprocal of
 792 the simplest rational between d'%r' and d%r.
 793
 794 \begin{code}
 795 approxRational x eps    =  simplest (x-eps) (x+eps)
 796         where simplest x y | y < x      =  simplest y x
 797                            | x == y     =  xr
 798                            | x > 0      =  simplest' n d n' d'
 799                            | y < 0      =  - simplest' (-n') d' (-n) d
 800                            | otherwise  =  0 :% 1
 801                                         where xr  = toRational x
 802                                               n   = numerator xr
 803                                               d   = denominator xr
 804                                               nd' = toRational y
 805                                               n'  = numerator nd'
 806                                               d'  = denominator nd'
 807
 808               simplest' n d n' d'       -- assumes 0 < n%d < n'%d'
 809                         | r == 0     =  q :% 1
 810                         | q /= q'    =  (q+1) :% 1
 811                         | otherwise  =  (q*n''+d'') :% n''
 812                                      where (q,r)      =  quotRem n d
 813                                            (q',r')    =  quotRem n' d'
 814                                            nd''       =  simplest' d' r' d r
 815                                            n''        =  numerator nd''
 816                                            d''        =  denominator nd''
 817 \end{code}
 818
 819
 820 \begin{code}
 821 instance  (Integral a)  => Ord (Ratio a)  where
 822     (x:%y) <= (x':%y')  =  x * y' <= x' * y
 823     (x:%y) <  (x':%y')  =  x * y' <  x' * y
 824
 825 instance  (Integral a)  => Num (Ratio a)  where
 826     (x:%y) + (x':%y')   =  reduce (x*y' + x'*y) (y*y')
 827     (x:%y) - (x':%y')   =  reduce (x*y' - x'*y) (y*y')
 828     (x:%y) * (x':%y')   =  reduce (x * x') (y * y')
 829     negate (x:%y)       =  (-x) :% y
 830     abs (x:%y)          =  abs x :% y
 831     signum (x:%y)       =  signum x :% 1
 832     fromInteger x       =  fromInteger x :% 1
 833
 834 instance  (Integral a)  => Real (Ratio a)  where
 835     toRational (x:%y)   =  toInteger x :% toInteger y
 836
 837 instance  (Integral a)  => Fractional (Ratio a)  where
 838     (x:%y) / (x':%y')   =  (x*y') % (y*x')
 839     recip (x:%y)        =  if x < 0 then (-y) :% (-x) else y :% x
 840     fromRational (x:%y) =  fromInteger x :% fromInteger y
 841
 842 instance  (Integral a)  => RealFrac (Ratio a)  where
 843     properFraction (x:%y) = (fromIntegral q, r:%y)
 844                             where (q,r) = quotRem x y
 845
 846 instance  (Integral a)  => Enum (Ratio a)  where
 847     enumFrom            =  iterate ((+)1)
 848     enumFromThen n m    =  iterate ((+)(m-n)) n
 849     toEnum n            =  fromIntegral n :% 1
 850     fromEnum            =  fromInteger . truncate
 851
 852 ratio_prec :: Int
 853 ratio_prec = 7
 854
 855 instance  (Integral a)  => Show (Ratio a)  where
 856     showsPrec p (x:%y)  =  showParen (p > ratio_prec)
 857                                (shows x . showString " % " . shows y)
 858 \end{code}
 859
 860 \begin{code}
 861 --Exported from std library Numeric, defined here to
 862 --avoid mut. rec. between PrelNum and Numeric.
 863 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
 864 showSigned showPos p x = if x < 0 then showParen (p > 6)
 865                                                  (showChar '-' . showPos (-x))
 866                                   else showPos x
 867
 868 showSignedInteger :: Int -> Integer -> ShowS
 869 showSignedInteger p n r
 870   = -- from HBC version; support code follows
 871     if n < 0 && p > 6 then '(':jtos n++(')':r) else jtos n ++ r
 872
 873 jtos :: Integer -> String
 874 jtos n
 875   = if n < 0 then
 876         '-' : jtos' (-n) []
 877     else
 878         jtos' n []
 879
 880 jtos' :: Integer -> String -> String
 881 jtos' n cs
 882   = if n < 10 then
 883         chr (fromInteger (n + ord_0)) : cs
 884     else
 885         jtos' q (chr (toInt r + (ord_0::Int)) : cs)
 886   where
 887     (q,r) = n `quotRem` 10
 888
 889 showFloat x  =  showString (formatRealFloat FFGeneric Nothing x)
 890
 891 -- These are the format types.  This type is not exported.
 892
 893 data FFFormat = FFExponent | FFFixed | FFGeneric --no need: deriving (Eq, Ord, Show)
 894
 895 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
 896 formatRealFloat fmt decs x = s
 897  where
 898   base = 10
 899   s = if isNaN x
 900       then "NaN"
 901       else
 902        if isInfinite x then
 903           if x < 0 then "-Infinity" else "Infinity"
 904        else
 905           if x < 0 || isNegativeZero x then
 906             '-':doFmt fmt (floatToDigits (toInteger base) (-x))
 907           else
 908             doFmt fmt (floatToDigits (toInteger base) x)
 909
 910   doFmt fmt (is, e) =
 911     let ds = map intToDigit is in
 912     case fmt of
 913      FFGeneric ->
 914       doFmt (if e <0 || e > 7 then FFExponent else FFFixed)
 915             (is,e)
 916      FFExponent ->
 917       case decs of
 918        Nothing ->
 919         let e' = if e==0 then 0 else e-1 in
 920         (case ds of
 921           [d]    -> d : ".0e"
 922           (d:ds) -> d : '.' : ds ++ "e") ++ show e'
 923        Just dec ->
 924         let dec' = max dec 1 in
 925         case is of
 926          [0] -> '0':'.':take dec' (repeat '0') ++ "e0"
 927          _ ->
 928           let
 929            (ei,is') = roundTo base (dec'+1) is
 930            d:ds = map intToDigit (if ei > 0 then init is' else is')
 931           in
 932           d:'.':ds ++ 'e':show (e-1+ei)
 933      FFFixed ->
 934       let
 935        mk0 ls = case ls of { "" -> "0" ; _ -> ls}
 936       in
 937       case decs of
 938        Nothing ->
 939          let
 940           f 0 s ds = mk0 (reverse s) ++ '.':mk0 ds
 941           f n s "" = f (n-1) ('0':s) ""
 942           f n s (d:ds) = f (n-1) (d:s) ds
 943          in
 944          f e "" ds
 945        Just dec ->
 946         let dec' = max dec 1 in
 947         if e >= 0 then
 948          let
 949           (ei,is') = roundTo base (dec' + e) is
 950           (ls,rs)  = splitAt (e+ei) (map intToDigit is')
 951          in
 952          mk0 ls ++ (if null rs then "" else '.':rs)
 953         else
 954          let
 955           (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
 956           d:ds = map intToDigit (if ei > 0 then is' else 0:is')
 957          in
 958          d : '.' : ds
 959
 960
 961 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
 962 roundTo base d is =
 963  let
 964   v = f d is
 965  in
 966  case v of
 967   (0,is) -> v
 968   (1,is) -> (1, 1:is)
 969  where
 970   b2 = base `div` 2
 971
 972   f n [] = (0, replicate n 0)
 973   f 0 (i:_) = (if i>=b2 then 1 else 0, [])
 974   f d (i:is) =
 975     let
 976      (c,ds) = f (d-1) is
 977      i' = c + i
 978     in
 979     if i' == base then (1,0:ds) else (0,i':ds)
 980
 981 --
 982 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
 983 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
 984 -- This version uses a much slower logarithm estimator. It should be improved.
 985
 986 -- This function returns a list of digits (Ints in [0..base-1]) and an
 987 -- exponent.
 988 --floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
 989 floatToDigits _ 0 = ([0], 0)
 990 floatToDigits base x =
 991  let
 992   (f0, e0) = decodeFloat x
 993   (minExp0, _) = floatRange x
 994   p = floatDigits x
 995   b = floatRadix x
 996   minExp = minExp0 - p -- the real minimum exponent
 997   -- Haskell requires that f be adjusted so denormalized numbers
 998   -- will have an impossibly low exponent.  Adjust for this.
 999   (f, e) =
1000    let n = minExp - e0 in
1001    if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
1002   (r, s, mUp, mDn) =
1003    if e >= 0 then
1004     let be = b^ e in
1005     if f == b^(p-1) then
1006       (f*be*b*2, 2*b, be*b, b)
1007     else
1008       (f*be*2, 2, be, be)
1009    else
1010     if e > minExp && f == b^(p-1) then
1011       (f*b*2, b^(-e+1)*2, b, 1)
1012     else
1013       (f*2, b^(-e)*2, 1, 1)
1014   k =
1015    let
1016     k0 =
1017      if b == 2 && base == 10 then
1018         -- logBase 10 2 is slightly bigger than 3/10 so
1019         -- the following will err on the low side.  Ignoring
1020         -- the fraction will make it err even more.
1021         -- Haskell promises that p-1 <= logBase b f < p.
1022         (p - 1 + e0) * 3 `div` 10
1023      else
1024         ceiling ((log (fromInteger (f+1)) + fromInt e * log (fromInteger b)) /
1025                   log (fromInteger base))
1026
1027     fixup n =
1028       if n >= 0 then
1029         if r + mUp <= expt base n * s then n else fixup (n+1)
1030       else
1031         if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
1032    in
1033    fixup k0
1034
1035   gen ds rn sN mUpN mDnN =
1036    let
1037     (dn, rn') = (rn * base) `divMod` sN
1038     mUpN' = mUpN * base
1039     mDnN' = mDnN * base
1040    in
1041    case (rn' < mDnN', rn' + mUpN' > sN) of
1042     (True,  False) -> dn : ds
1043     (False, True)  -> dn+1 : ds
1044     (True,  True)  -> if rn' * 2 < sN then dn : ds else dn+1 : ds
1045     (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
1046
1047   rds =
1048    if k >= 0 then
1049       gen [] r (s * expt base k) mUp mDn
1050    else
1051      let bk = expt base (-k) in
1052      gen [] (r * bk) s (mUp * bk) (mDn * bk)
1053  in
1054  (map toInt (reverse rds), k)
1055
1056 \end{code}
1057
1058 @showRational@ converts a Rational to a string that looks like a
1059 floating point number, but without converting to any floating type
1060 (because of the possible overflow).
1061
1062 From/by Lennart, 94/09/26
1063
1064 \begin{code}
1065 showRational :: Int -> Rational -> String
1066 showRational n r =
1067     if r == 0 then
1068         "0.0"
1069     else
1070         let (r', e) = normalize r
1071         in  prR n r' e
1072
1073 startExpExp = 4 :: Int
1074
1075 -- make sure 1 <= r < 10
1076 normalize :: Rational -> (Rational, Int)
1077 normalize r = if r < 1 then
1078                   case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
1079               else
1080                   norm startExpExp r 0
1081         where norm :: Int -> Rational -> Int -> (Rational, Int)
1082               -- Invariant: r*10^e == original r
1083               norm 0  r e = (r, e)
1084               norm ee r e =
1085                 let n = 10^ee
1086                     tn = 10^n
1087                 in  if r >= tn then norm ee (r/tn) (e+n) else norm (ee-1) r e
1088
1089 drop0 "" = ""
1090 drop0 (c:cs) = c : reverse (dropWhile (=='0') (reverse cs))
1091
1092 prR :: Int -> Rational -> Int -> String
1093 prR n r e | r <  1  = prR n (r*10) (e-1)                -- final adjustment
1094 prR n r e | r >= 10 = prR n (r/10) (e+1)
1095 prR n r e0 =
1096         let s = show ((round (r * 10^n))::Integer)
1097             e = e0+1
1098         in  if e > 0 && e < 8 then
1099                 take e s ++ "." ++ drop0 (drop e s)
1100             else if e <= 0 && e > -3 then
1101                 "0." ++ take (-e) (repeat '0') ++ drop0 s
1102             else
1103                 head s : "."++ drop0 (tail s) ++ "e" ++ show e0
1104 \end{code}
1105
1106
1107 [In response to a request for documentation of how fromRational works,
1108 Joe Fasel writes:] A quite reasonable request!  This code was added to
1109 the Prelude just before the 1.2 release, when Lennart, working with an
1110 early version of hbi, noticed that (read . show) was not the identity
1111 for floating-point numbers.  (There was a one-bit error about half the
1112 time.)  The original version of the conversion function was in fact
1113 simply a floating-point divide, as you suggest above. The new version
1114 is, I grant you, somewhat denser.
1115
1116 Unfortunately, Joe's code doesn't work!  Here's an example:
1117
1118 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
1119
1120 This program prints
1121         0.0000000000000000
1122 instead of
1123         1.8217369128763981e-300
1124
1125 Lennart's code follows, and it works...
1126
1127 \begin{pseudocode}
1128 {-# GENERATE_SPECS fromRational__ a{Double#,Double} #-}
1129 fromRat :: (RealFloat a) => Rational -> a
1130 fromRat x = x'
1131         where x' = f e
1132
1133 --              If the exponent of the nearest floating-point number to x
1134 --              is e, then the significand is the integer nearest xb^(-e),
1135 --              where b is the floating-point radix.  We start with a good
1136 --              guess for e, and if it is correct, the exponent of the
1137 --              floating-point number we construct will again be e.  If
1138 --              not, one more iteration is needed.
1139
1140               f e   = if e' == e then y else f e'
1141                       where y      = encodeFloat (round (x * (1 % b)^^e)) e
1142                             (_,e') = decodeFloat y
1143               b     = floatRadix x'
1144
1145 --              We obtain a trial exponent by doing a floating-point
1146 --              division of x's numerator by its denominator.  The
1147 --              result of this division may not itself be the ultimate
1148 --              result, because of an accumulation of three rounding
1149 --              errors.
1150
1151               (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
1152                                         / fromInteger (denominator x))
1153 \end{pseudocode}
1154
1155 Now, here's Lennart's code.
1156
1157 \begin{code}
1158 {-# SPECIALISE fromRat ::
1159         Rational -> Double,
1160         Rational -> Float #-}
1161
1162 --fromRat :: (RealFloat a) => Rational -> a
1163 fromRat x =
1164     if x == 0 then encodeFloat 0 0              -- Handle exceptional cases
1165     else if x < 0 then - fromRat' (-x)          -- first.
1166     else fromRat' x
1167
1168 -- Conversion process:
1169 -- Scale the rational number by the RealFloat base until
1170 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
1171 -- Then round the rational to an Integer and encode it with the exponent
1172 -- that we got from the scaling.
1173 -- To speed up the scaling process we compute the log2 of the number to get
1174 -- a first guess of the exponent.
1175
1176 fromRat' :: (RealFloat a) => Rational -> a
1177 fromRat' x = r
1178   where b = floatRadix r
1179         p = floatDigits r
1180         (minExp0, _) = floatRange r
1181         minExp = minExp0 - p            -- the real minimum exponent
1182         xMin = toRational (expt b (p-1))
1183         xMax = toRational (expt b p)
1184         p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
1185         f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
1186         (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
1187         r = encodeFloat (round x') p'
1188
1189 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
1190 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
1191 scaleRat b minExp xMin xMax p x
1192     | p <= minExp = (x, p)
1193     | x >= xMax   = scaleRat b minExp xMin xMax (p+1) (x/b)
1194     | x < xMin    = scaleRat b minExp xMin xMax (p-1) (x*b)
1195     | otherwise   = (x, p)
1196
1197 -- Exponentiation with a cache for the most common numbers.
1198 minExpt = 0::Int
1199 maxExpt = 1100::Int
1200 expt :: Integer -> Int -> Integer
1201 expt base n =
1202     if base == 2 && n >= minExpt && n <= maxExpt then
1203         expts!n
1204     else
1205         base^n
1206 expts :: Array Int Integer
1207 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
1208
1209 -- Compute the (floor of the) log of i in base b.
1210 -- Simplest way would be just divide i by b until it's smaller then b, but that would
1211 -- be very slow!  We are just slightly more clever.
1212 integerLogBase :: Integer -> Integer -> Int
1213 integerLogBase b i =
1214      if i < b then
1215         0
1216      else
1217         -- Try squaring the base first to cut down the number of divisions.
1218         let l = 2 * integerLogBase (b*b) i
1219             doDiv :: Integer -> Int -> Int
1220             doDiv i l = if i < b then l else doDiv (i `div` b) (l+1)
1221         in  doDiv (i `div` (b^l)) l
1222 \end{code}
1223
1224 %*********************************************************
1225 %*                                                      *
1226 \subsection{Numeric primops}
1227 %*                                                      *
1228 %*********************************************************
1229
1230 Definitions of the boxed PrimOps; these will be
1231 used in the case of partial applications, etc.
1232
1233 \begin{code}
1234 plusFloat   (F# x) (F# y) = F# (plusFloat# x y)
1235 minusFloat  (F# x) (F# y) = F# (minusFloat# x y)
1236 timesFloat  (F# x) (F# y) = F# (timesFloat# x y)
1237 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
1238 negateFloat (F# x)        = F# (negateFloat# x)
1239
1240 gtFloat     (F# x) (F# y) = gtFloat# x y
1241 geFloat     (F# x) (F# y) = geFloat# x y
1242 eqFloat     (F# x) (F# y) = eqFloat# x y
1243 neFloat     (F# x) (F# y) = neFloat# x y
1244 ltFloat     (F# x) (F# y) = ltFloat# x y
1245 leFloat     (F# x) (F# y) = leFloat# x y
1246
1247 float2Int   (F# x) = I# (float2Int# x)
1248 int2Float   (I# x) = F# (int2Float# x)
1249
1250 expFloat    (F# x) = F# (expFloat# x)
1251 logFloat    (F# x) = F# (logFloat# x)
1252 sqrtFloat   (F# x) = F# (sqrtFloat# x)
1253 sinFloat    (F# x) = F# (sinFloat# x)
1254 cosFloat    (F# x) = F# (cosFloat# x)
1255 tanFloat    (F# x) = F# (tanFloat# x)
1256 asinFloat   (F# x) = F# (asinFloat# x)
1257 acosFloat   (F# x) = F# (acosFloat# x)
1258 atanFloat   (F# x) = F# (atanFloat# x)
1259 sinhFloat   (F# x) = F# (sinhFloat# x)
1260 coshFloat   (F# x) = F# (coshFloat# x)
1261 tanhFloat   (F# x) = F# (tanhFloat# x)
1262
1263 powerFloat  (F# x) (F# y) = F# (powerFloat# x y)
1264
1265 -- definitions of the boxed PrimOps; these will be
1266 -- used in the case of partial applications, etc.
1267
1268 plusDouble   (D# x) (D# y) = D# (x +## y)
1269 minusDouble  (D# x) (D# y) = D# (x -## y)
1270 timesDouble  (D# x) (D# y) = D# (x *## y)
1271 divideDouble (D# x) (D# y) = D# (x /## y)
1272 negateDouble (D# x)        = D# (negateDouble# x)
1273
1274 gtDouble    (D# x) (D# y) = x >## y
1275 geDouble    (D# x) (D# y) = x >=## y
1276 eqDouble    (D# x) (D# y) = x ==## y
1277 neDouble    (D# x) (D# y) = x /=## y
1278 ltDouble    (D# x) (D# y) = x <## y
1279 leDouble    (D# x) (D# y) = x <=## y
1280
1281 double2Int   (D# x) = I# (double2Int#   x)
1282 int2Double   (I# x) = D# (int2Double#   x)
1283 double2Float (D# x) = F# (double2Float# x)
1284 float2Double (F# x) = D# (float2Double# x)
1285
1286 expDouble    (D# x) = D# (expDouble# x)
1287 logDouble    (D# x) = D# (logDouble# x)
1288 sqrtDouble   (D# x) = D# (sqrtDouble# x)
1289 sinDouble    (D# x) = D# (sinDouble# x)
1290 cosDouble    (D# x) = D# (cosDouble# x)
1291 tanDouble    (D# x) = D# (tanDouble# x)
1292 asinDouble   (D# x) = D# (asinDouble# x)
1293 acosDouble   (D# x) = D# (acosDouble# x)
1294 atanDouble   (D# x) = D# (atanDouble# x)
1295 sinhDouble   (D# x) = D# (sinhDouble# x)
1296 coshDouble   (D# x) = D# (coshDouble# x)
1297 tanhDouble   (D# x) = D# (tanhDouble# x)
1298
1299 powerDouble  (D# x) (D# y) = D# (x **## y)
1300 \end{code}