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