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