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)) +
1025                  fromInt e * log (fromInteger b)) /
1026                   fromInt e * log (fromInteger b))
1027
1028     fixup n =
1029       if n >= 0 then
1030         if r + mUp <= expt base n * s then n else fixup (n+1)
1031       else
1032         if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
1033    in
1034    fixup k0
1035
1036   gen ds rn sN mUpN mDnN =
1037    let
1038     (dn, rn') = (rn * base) `divMod` sN
1039     mUpN' = mUpN * base
1040     mDnN' = mDnN * base
1041    in
1042    case (rn' < mDnN', rn' + mUpN' > sN) of
1043     (True,  False) -> dn : ds
1044     (False, True)  -> dn+1 : ds
1045     (True,  True)  -> if rn' * 2 < sN then dn : ds else dn+1 : ds
1046     (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
1047
1048   rds =
1049    if k >= 0 then
1050       gen [] r (s * expt base k) mUp mDn
1051    else
1052      let bk = expt base (-k) in
1053      gen [] (r * bk) s (mUp * bk) (mDn * bk)
1054  in
1055  (map toInt (reverse rds), k)
1056
1057 \end{code}
1058
1059 @showRational@ converts a Rational to a string that looks like a
1060 floating point number, but without converting to any floating type
1061 (because of the possible overflow).
1062
1063 From/by Lennart, 94/09/26
1064
1065 \begin{code}
1066 showRational :: Int -> Rational -> String
1067 showRational n r =
1068     if r == 0 then
1069         "0.0"
1070     else
1071         let (r', e) = normalize r
1072         in  prR n r' e
1073
1074 startExpExp = 4 :: Int
1075
1076 -- make sure 1 <= r < 10
1077 normalize :: Rational -> (Rational, Int)
1078 normalize r = if r < 1 then
1079                   case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
1080               else
1081                   norm startExpExp r 0
1082         where norm :: Int -> Rational -> Int -> (Rational, Int)
1083               -- Invariant: r*10^e == original r
1084               norm 0  r e = (r, e)
1085               norm ee r e =
1086                 let n = 10^ee
1087                     tn = 10^n
1088                 in  if r >= tn then norm ee (r/tn) (e+n) else norm (ee-1) r e
1089
1090 drop0 "" = ""
1091 drop0 (c:cs) = c : reverse (dropWhile (=='0') (reverse cs))
1092
1093 prR :: Int -> Rational -> Int -> String
1094 prR n r e | r <  1  = prR n (r*10) (e-1)                -- final adjustment
1095 prR n r e | r >= 10 = prR n (r/10) (e+1)
1096 prR n r e0 =
1097         let s = show ((round (r * 10^n))::Integer)
1098             e = e0+1
1099         in  if e > 0 && e < 8 then
1100                 take e s ++ "." ++ drop0 (drop e s)
1101             else if e <= 0 && e > -3 then
1102                 "0." ++ take (-e) (repeat '0') ++ drop0 s
1103             else
1104                 head s : "."++ drop0 (tail s) ++ "e" ++ show e0
1105 \end{code}
1106
1107
1108 [In response to a request for documentation of how fromRational works,
1109 Joe Fasel writes:] A quite reasonable request!  This code was added to
1110 the Prelude just before the 1.2 release, when Lennart, working with an
1111 early version of hbi, noticed that (read . show) was not the identity
1112 for floating-point numbers.  (There was a one-bit error about half the
1113 time.)  The original version of the conversion function was in fact
1114 simply a floating-point divide, as you suggest above. The new version
1115 is, I grant you, somewhat denser.
1116
1117 Unfortunately, Joe's code doesn't work!  Here's an example:
1118
1119 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
1120
1121 This program prints
1122         0.0000000000000000
1123 instead of
1124         1.8217369128763981e-300
1125
1126 Lennart's code follows, and it works...
1127
1128 \begin{pseudocode}
1129 {-# GENERATE_SPECS fromRational__ a{Double#,Double} #-}
1130 fromRat :: (RealFloat a) => Rational -> a
1131 fromRat x = x'
1132         where x' = f e
1133
1134 --              If the exponent of the nearest floating-point number to x
1135 --              is e, then the significand is the integer nearest xb^(-e),
1136 --              where b is the floating-point radix.  We start with a good
1137 --              guess for e, and if it is correct, the exponent of the
1138 --              floating-point number we construct will again be e.  If
1139 --              not, one more iteration is needed.
1140
1141               f e   = if e' == e then y else f e'
1142                       where y      = encodeFloat (round (x * (1 % b)^^e)) e
1143                             (_,e') = decodeFloat y
1144               b     = floatRadix x'
1145
1146 --              We obtain a trial exponent by doing a floating-point
1147 --              division of x's numerator by its denominator.  The
1148 --              result of this division may not itself be the ultimate
1149 --              result, because of an accumulation of three rounding
1150 --              errors.
1151
1152               (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
1153                                         / fromInteger (denominator x))
1154 \end{pseudocode}
1155
1156 Now, here's Lennart's code.
1157
1158 \begin{code}
1159 {-# SPECIALISE fromRat ::
1160         Rational -> Double,
1161         Rational -> Float #-}
1162
1163 --fromRat :: (RealFloat a) => Rational -> a
1164 fromRat x =
1165     if x == 0 then encodeFloat 0 0              -- Handle exceptional cases
1166     else if x < 0 then - fromRat' (-x)          -- first.
1167     else fromRat' x
1168
1169 -- Conversion process:
1170 -- Scale the rational number by the RealFloat base until
1171 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
1172 -- Then round the rational to an Integer and encode it with the exponent
1173 -- that we got from the scaling.
1174 -- To speed up the scaling process we compute the log2 of the number to get
1175 -- a first guess of the exponent.
1176
1177 fromRat' :: (RealFloat a) => Rational -> a
1178 fromRat' x = r
1179   where b = floatRadix r
1180         p = floatDigits r
1181         (minExp0, _) = floatRange r
1182         minExp = minExp0 - p            -- the real minimum exponent
1183         xMin = toRational (expt b (p-1))
1184         xMax = toRational (expt b p)
1185         p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
1186         f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
1187         (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
1188         r = encodeFloat (round x') p'
1189
1190 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
1191 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
1192 scaleRat b minExp xMin xMax p x =
1193     if p <= minExp then
1194         (x, p)
1195     else if x >= xMax then
1196         scaleRat b minExp xMin xMax (p+1) (x/b)
1197     else if x < xMin  then
1198         scaleRat b minExp xMin xMax (p-1) (x*b)
1199     else
1200         (x, p)
1201
1202 -- Exponentiation with a cache for the most common numbers.
1203 minExpt = 0::Int
1204 maxExpt = 1100::Int
1205 expt :: Integer -> Int -> Integer
1206 expt base n =
1207     if base == 2 && n >= minExpt && n <= maxExpt then
1208         expts!n
1209     else
1210         base^n
1211 expts :: Array Int Integer
1212 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
1213
1214 -- Compute the (floor of the) log of i in base b.
1215 -- Simplest way would be just divide i by b until it's smaller then b, but that would
1216 -- be very slow!  We are just slightly more clever.
1217 integerLogBase :: Integer -> Integer -> Int
1218 integerLogBase b i =
1219      if i < b then
1220         0
1221      else
1222         -- Try squaring the base first to cut down the number of divisions.
1223         let l = 2 * integerLogBase (b*b) i
1224             doDiv :: Integer -> Int -> Int
1225             doDiv i l = if i < b then l else doDiv (i `div` b) (l+1)
1226         in  doDiv (i `div` (b^l)) l
1227 \end{code}
1228
1229 %*********************************************************
1230 %*                                                      *
1231 \subsection{Numeric primops}
1232 %*                                                      *
1233 %*********************************************************
1234
1235 Definitions of the boxed PrimOps; these will be
1236 used in the case of partial applications, etc.
1237
1238 \begin{code}
1239 plusFloat   (F# x) (F# y) = F# (plusFloat# x y)
1240 minusFloat  (F# x) (F# y) = F# (minusFloat# x y)
1241 timesFloat  (F# x) (F# y) = F# (timesFloat# x y)
1242 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
1243 negateFloat (F# x)        = F# (negateFloat# x)
1244
1245 gtFloat     (F# x) (F# y) = gtFloat# x y
1246 geFloat     (F# x) (F# y) = geFloat# x y
1247 eqFloat     (F# x) (F# y) = eqFloat# x y
1248 neFloat     (F# x) (F# y) = neFloat# x y
1249 ltFloat     (F# x) (F# y) = ltFloat# x y
1250 leFloat     (F# x) (F# y) = leFloat# x y
1251
1252 float2Int   (F# x) = I# (float2Int# x)
1253 int2Float   (I# x) = F# (int2Float# x)
1254
1255 expFloat    (F# x) = F# (expFloat# x)
1256 logFloat    (F# x) = F# (logFloat# x)
1257 sqrtFloat   (F# x) = F# (sqrtFloat# x)
1258 sinFloat    (F# x) = F# (sinFloat# x)
1259 cosFloat    (F# x) = F# (cosFloat# x)
1260 tanFloat    (F# x) = F# (tanFloat# x)
1261 asinFloat   (F# x) = F# (asinFloat# x)
1262 acosFloat   (F# x) = F# (acosFloat# x)
1263 atanFloat   (F# x) = F# (atanFloat# x)
1264 sinhFloat   (F# x) = F# (sinhFloat# x)
1265 coshFloat   (F# x) = F# (coshFloat# x)
1266 tanhFloat   (F# x) = F# (tanhFloat# x)
1267
1268 powerFloat  (F# x) (F# y) = F# (powerFloat# x y)
1269
1270 -- definitions of the boxed PrimOps; these will be
1271 -- used in the case of partial applications, etc.
1272
1273 plusDouble   (D# x) (D# y) = D# (x +## y)
1274 minusDouble  (D# x) (D# y) = D# (x -## y)
1275 timesDouble  (D# x) (D# y) = D# (x *## y)
1276 divideDouble (D# x) (D# y) = D# (x /## y)
1277 negateDouble (D# x)        = D# (negateDouble# x)
1278
1279 gtDouble    (D# x) (D# y) = x >## y
1280 geDouble    (D# x) (D# y) = x >=## y
1281 eqDouble    (D# x) (D# y) = x ==## y
1282 neDouble    (D# x) (D# y) = x /=## y
1283 ltDouble    (D# x) (D# y) = x <## y
1284 leDouble    (D# x) (D# y) = x <=## y
1285
1286 double2Int   (D# x) = I# (double2Int#   x)
1287 int2Double   (I# x) = D# (int2Double#   x)
1288 double2Float (D# x) = F# (double2Float# x)
1289 float2Double (F# x) = D# (float2Double# x)
1290
1291 expDouble    (D# x) = D# (expDouble# x)
1292 logDouble    (D# x) = D# (logDouble# x)
1293 sqrtDouble   (D# x) = D# (sqrtDouble# x)
1294 sinDouble    (D# x) = D# (sinDouble# x)
1295 cosDouble    (D# x) = D# (cosDouble# x)
1296 tanDouble    (D# x) = D# (tanDouble# x)
1297 asinDouble   (D# x) = D# (asinDouble# x)
1298 acosDouble   (D# x) = D# (acosDouble# x)
1299 atanDouble   (D# x) = D# (atanDouble# x)
1300 sinhDouble   (D# x) = D# (sinhDouble# x)
1301 coshDouble   (D# x) = D# (coshDouble# x)
1302 tanhDouble   (D# x) = D# (tanhDouble# x)
1303
1304 powerDouble  (D# x) (D# y) = D# (x **## y)
1305 \end{code}