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' q (chr (toInt r + (ord_0::Int)) : cs)
 852   where
 853     (q,r) = n `quotRem` 10
 854
 855 showFloat x  =  showString (formatRealFloat FFGeneric Nothing x)
 856
 857 -- These are the format types.  This type is not exported.
 858
 859 data FFFormat = FFExponent | FFFixed | FFGeneric --no need: deriving (Eq, Ord, Show)
 860
 861 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
 862 formatRealFloat fmt decs x = s
 863  where
 864   base = 10
 865   s = if isNaN x
 866       then "NaN"
 867       else
 868        if isInfinite x then
 869           if x < 0 then "-Infinity" else "Infinity"
 870        else
 871           if x < 0 || isNegativeZero x then
 872             '-':doFmt fmt (floatToDigits (toInteger base) (-x))
 873           else
 874             doFmt fmt (floatToDigits (toInteger base) x)
 875
 876   doFmt fmt (is, e) =
 877     let ds = map intToDigit is in
 878     case fmt of
 879      FFGeneric ->
 880       doFmt (if e <0 || e > 7 then FFExponent else FFFixed)
 881             (is,e)
 882      FFExponent ->
 883       case decs of
 884        Nothing ->
 885         let e' = if e==0 then 0 else e-1 in
 886         (case ds of
 887           [d]    -> d : ".0e"
 888           (d:ds) -> d : '.' : ds ++ "e") ++ show e'
 889        Just dec ->
 890         let dec' = max dec 1 in
 891         case is of
 892          [0] -> '0':'.':take dec' (repeat '0') ++ "e0"
 893          _ ->
 894           let
 895            (ei,is') = roundTo base (dec'+1) is
 896            d:ds = map intToDigit (if ei > 0 then init is' else is')
 897           in
 898           d:'.':ds ++ 'e':show (e-1+ei)
 899      FFFixed ->
 900       let
 901        mk0 ls = case ls of { "" -> "0" ; _ -> ls}
 902       in
 903       case decs of
 904        Nothing ->
 905          let
 906           f 0 s ds = mk0 (reverse s) ++ '.':mk0 ds
 907           f n s "" = f (n-1) ('0':s) ""
 908           f n s (d:ds) = f (n-1) (d:s) ds
 909          in
 910          f e "" ds
 911        Just dec ->
 912         let dec' = max dec 1 in
 913         if e >= 0 then
 914          let
 915           (ei,is') = roundTo base (dec' + e) is
 916           (ls,rs)  = splitAt (e+ei) (map intToDigit is')
 917          in
 918          mk0 ls ++ (if null rs then "" else '.':rs)
 919         else
 920          let
 921           (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
 922           d:ds = map intToDigit (if ei > 0 then is' else 0:is')
 923          in
 924          d : '.' : ds
 925
 926
 927 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
 928 roundTo base d is =
 929  let
 930   v = f d is
 931  in
 932  case v of
 933   (0,is) -> v
 934   (1,is) -> (1, 1:is)
 935  where
 936   b2 = base `div` 2
 937
 938   f n [] = (0, replicate n 0)
 939   f 0 (i:_) = (if i>=b2 then 1 else 0, [])
 940   f d (i:is) =
 941     let
 942      (c,ds) = f (d-1) is
 943      i' = c + i
 944     in
 945     if i' == base then (1,0:ds) else (0,i':ds)
 946
 947 --
 948 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
 949 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
 950 -- This version uses a much slower logarithm estimator. It should be improved.
 951
 952 -- This function returns a list of digits (Ints in [0..base-1]) and an
 953 -- exponent.
 954 --floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
 955 floatToDigits _ 0 = ([0], 0)
 956 floatToDigits base x =
 957  let
 958   (f0, e0) = decodeFloat x
 959   (minExp0, _) = floatRange x
 960   p = floatDigits x
 961   b = floatRadix x
 962   minExp = minExp0 - p -- the real minimum exponent
 963   -- Haskell requires that f be adjusted so denormalized numbers
 964   -- will have an impossibly low exponent.  Adjust for this.
 965   (f, e) =
 966    let n = minExp - e0 in
 967    if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
 968   (r, s, mUp, mDn) =
 969    if e >= 0 then
 970     let be = b^ e in
 971     if f == b^(p-1) then
 972       (f*be*b*2, 2*b, be*b, b)
 973     else
 974       (f*be*2, 2, be, be)
 975    else
 976     if e > minExp && f == b^(p-1) then
 977       (f*b*2, b^(-e+1)*2, b, 1)
 978     else
 979       (f*2, b^(-e)*2, 1, 1)
 980   k =
 981    let
 982     k0 =
 983      if b == 2 && base == 10 then
 984         -- logBase 10 2 is slightly bigger than 3/10 so
 985         -- the following will err on the low side.  Ignoring
 986         -- the fraction will make it err even more.
 987         -- Haskell promises that p-1 <= logBase b f < p.
 988         (p - 1 + e0) * 3 `div` 10
 989      else
 990         ceiling ((log (fromInteger (f+1)) +
 991                  fromInt e * log (fromInteger b)) /
 992                   fromInt e * log (fromInteger b))
 993
 994     fixup n =
 995       if n >= 0 then
 996         if r + mUp <= expt base n * s then n else fixup (n+1)
 997       else
 998         if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
 999    in
1000    fixup k0
1001
1002   gen ds rn sN mUpN mDnN =
1003    let
1004     (dn, rn') = (rn * base) `divMod` sN
1005     mUpN' = mUpN * base
1006     mDnN' = mDnN * base
1007    in
1008    case (rn' < mDnN', rn' + mUpN' > sN) of
1009     (True,  False) -> dn : ds
1010     (False, True)  -> dn+1 : ds
1011     (True,  True)  -> if rn' * 2 < sN then dn : ds else dn+1 : ds
1012     (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
1013
1014   rds =
1015    if k >= 0 then
1016       gen [] r (s * expt base k) mUp mDn
1017    else
1018      let bk = expt base (-k) in
1019      gen [] (r * bk) s (mUp * bk) (mDn * bk)
1020  in
1021  (map toInt (reverse rds), k)
1022
1023 \end{code}
1024
1025 @showRational@ converts a Rational to a string that looks like a
1026 floating point number, but without converting to any floating type
1027 (because of the possible overflow).
1028
1029 From/by Lennart, 94/09/26
1030
1031 \begin{code}
1032 showRational :: Int -> Rational -> String
1033 showRational n r =
1034     if r == 0 then
1035         "0.0"
1036     else
1037         let (r', e) = normalize r
1038         in  prR n r' e
1039
1040 startExpExp = 4 :: Int
1041
1042 -- make sure 1 <= r < 10
1043 normalize :: Rational -> (Rational, Int)
1044 normalize r = if r < 1 then
1045                   case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
1046               else
1047                   norm startExpExp r 0
1048         where norm :: Int -> Rational -> Int -> (Rational, Int)
1049               -- Invariant: r*10^e == original r
1050               norm 0  r e = (r, e)
1051               norm ee r e =
1052                 let n = 10^ee
1053                     tn = 10^n
1054                 in  if r >= tn then norm ee (r/tn) (e+n) else norm (ee-1) r e
1055
1056 drop0 "" = ""
1057 drop0 (c:cs) = c : reverse (dropWhile (=='0') (reverse cs))
1058
1059 prR :: Int -> Rational -> Int -> String
1060 prR n r e | r <  1  = prR n (r*10) (e-1)                -- final adjustment
1061 prR n r e | r >= 10 = prR n (r/10) (e+1)
1062 prR n r e0 =
1063         let s = show ((round (r * 10^n))::Integer)
1064             e = e0+1
1065         in  if e > 0 && e < 8 then
1066                 take e s ++ "." ++ drop0 (drop e s)
1067             else if e <= 0 && e > -3 then
1068                 "0." ++ take (-e) (repeat '0') ++ drop0 s
1069             else
1070                 head s : "."++ drop0 (tail s) ++ "e" ++ show e0
1071 \end{code}
1072
1073
1074 [In response to a request for documentation of how fromRational works,
1075 Joe Fasel writes:] A quite reasonable request!  This code was added to
1076 the Prelude just before the 1.2 release, when Lennart, working with an
1077 early version of hbi, noticed that (read . show) was not the identity
1078 for floating-point numbers.  (There was a one-bit error about half the
1079 time.)  The original version of the conversion function was in fact
1080 simply a floating-point divide, as you suggest above. The new version
1081 is, I grant you, somewhat denser.
1082
1083 Unfortunately, Joe's code doesn't work!  Here's an example:
1084
1085 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
1086
1087 This program prints
1088         0.0000000000000000
1089 instead of
1090         1.8217369128763981e-300
1091
1092 Lennart's code follows, and it works...
1093
1094 \begin{pseudocode}
1095 {-# GENERATE_SPECS fromRational__ a{Double#,Double} #-}
1096 fromRat :: (RealFloat a) => Rational -> a
1097 fromRat x = x'
1098         where x' = f e
1099
1100 --              If the exponent of the nearest floating-point number to x
1101 --              is e, then the significand is the integer nearest xb^(-e),
1102 --              where b is the floating-point radix.  We start with a good
1103 --              guess for e, and if it is correct, the exponent of the
1104 --              floating-point number we construct will again be e.  If
1105 --              not, one more iteration is needed.
1106
1107               f e   = if e' == e then y else f e'
1108                       where y      = encodeFloat (round (x * (1 % b)^^e)) e
1109                             (_,e') = decodeFloat y
1110               b     = floatRadix x'
1111
1112 --              We obtain a trial exponent by doing a floating-point
1113 --              division of x's numerator by its denominator.  The
1114 --              result of this division may not itself be the ultimate
1115 --              result, because of an accumulation of three rounding
1116 --              errors.
1117
1118               (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
1119                                         / fromInteger (denominator x))
1120 \end{pseudocode}
1121
1122 Now, here's Lennart's code.
1123
1124 \begin{code}
1125 --fromRat :: (RealFloat a) => Rational -> a
1126 fromRat x =
1127     if x == 0 then encodeFloat 0 0              -- Handle exceptional cases
1128     else if x < 0 then - fromRat' (-x)          -- first.
1129     else fromRat' x
1130
1131 -- Conversion process:
1132 -- Scale the rational number by the RealFloat base until
1133 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
1134 -- Then round the rational to an Integer and encode it with the exponent
1135 -- that we got from the scaling.
1136 -- To speed up the scaling process we compute the log2 of the number to get
1137 -- a first guess of the exponent.
1138
1139 fromRat' :: (RealFloat a) => Rational -> a
1140 fromRat' x = r
1141   where b = floatRadix r
1142         p = floatDigits r
1143         (minExp0, _) = floatRange r
1144         minExp = minExp0 - p            -- the real minimum exponent
1145         xMin = toRational (expt b (p-1))
1146         xMax = toRational (expt b p)
1147         p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
1148         f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
1149         (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
1150         r = encodeFloat (round x') p'
1151
1152 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
1153 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
1154 scaleRat b minExp xMin xMax p x =
1155     if p <= minExp then
1156         (x, p)
1157     else if x >= xMax then
1158         scaleRat b minExp xMin xMax (p+1) (x/b)
1159     else if x < xMin  then
1160         scaleRat b minExp xMin xMax (p-1) (x*b)
1161     else
1162         (x, p)
1163
1164 -- Exponentiation with a cache for the most common numbers.
1165 minExpt = 0::Int
1166 maxExpt = 1100::Int
1167 expt :: Integer -> Int -> Integer
1168 expt base n =
1169     if base == 2 && n >= minExpt && n <= maxExpt then
1170         expts!n
1171     else
1172         base^n
1173 expts :: Array Int Integer
1174 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
1175
1176 -- Compute the (floor of the) log of i in base b.
1177 -- Simplest way would be just divide i by b until it's smaller then b, but that would
1178 -- be very slow!  We are just slightly more clever.
1179 integerLogBase :: Integer -> Integer -> Int
1180 integerLogBase b i =
1181      if i < b then
1182         0
1183      else
1184         -- Try squaring the base first to cut down the number of divisions.
1185         let l = 2 * integerLogBase (b*b) i
1186             doDiv :: Integer -> Int -> Int
1187             doDiv i l = if i < b then l else doDiv (i `div` b) (l+1)
1188         in  doDiv (i `div` (b^l)) l
1189 \end{code}
1190
1191 %*********************************************************
1192 %*                                                      *
1193 \subsection{Numeric primops}
1194 %*                                                      *
1195 %*********************************************************
1196
1197 Definitions of the boxed PrimOps; these will be
1198 used in the case of partial applications, etc.
1199
1200 \begin{code}
1201 plusFloat   (F# x) (F# y) = F# (plusFloat# x y)
1202 minusFloat  (F# x) (F# y) = F# (minusFloat# x y)
1203 timesFloat  (F# x) (F# y) = F# (timesFloat# x y)
1204 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
1205 negateFloat (F# x)        = F# (negateFloat# x)
1206
1207 gtFloat     (F# x) (F# y) = gtFloat# x y
1208 geFloat     (F# x) (F# y) = geFloat# x y
1209 eqFloat     (F# x) (F# y) = eqFloat# x y
1210 neFloat     (F# x) (F# y) = neFloat# x y
1211 ltFloat     (F# x) (F# y) = ltFloat# x y
1212 leFloat     (F# x) (F# y) = leFloat# x y
1213
1214 float2Int   (F# x) = I# (float2Int# x)
1215 int2Float   (I# x) = F# (int2Float# x)
1216
1217 expFloat    (F# x) = F# (expFloat# x)
1218 logFloat    (F# x) = F# (logFloat# x)
1219 sqrtFloat   (F# x) = F# (sqrtFloat# x)
1220 sinFloat    (F# x) = F# (sinFloat# x)
1221 cosFloat    (F# x) = F# (cosFloat# x)
1222 tanFloat    (F# x) = F# (tanFloat# x)
1223 asinFloat   (F# x) = F# (asinFloat# x)
1224 acosFloat   (F# x) = F# (acosFloat# x)
1225 atanFloat   (F# x) = F# (atanFloat# x)
1226 sinhFloat   (F# x) = F# (sinhFloat# x)
1227 coshFloat   (F# x) = F# (coshFloat# x)
1228 tanhFloat   (F# x) = F# (tanhFloat# x)
1229
1230 powerFloat  (F# x) (F# y) = F# (powerFloat# x y)
1231
1232 -- definitions of the boxed PrimOps; these will be
1233 -- used in the case of partial applications, etc.
1234
1235 plusDouble   (D# x) (D# y) = D# (x +## y)
1236 minusDouble  (D# x) (D# y) = D# (x -## y)
1237 timesDouble  (D# x) (D# y) = D# (x *## y)
1238 divideDouble (D# x) (D# y) = D# (x /## y)
1239 negateDouble (D# x)        = D# (negateDouble# x)
1240
1241 gtDouble    (D# x) (D# y) = x >## y
1242 geDouble    (D# x) (D# y) = x >=## y
1243 eqDouble    (D# x) (D# y) = x ==## y
1244 neDouble    (D# x) (D# y) = x /=## y
1245 ltDouble    (D# x) (D# y) = x <## y
1246 leDouble    (D# x) (D# y) = x <=## y
1247
1248 double2Int   (D# x) = I# (double2Int#   x)
1249 int2Double   (I# x) = D# (int2Double#   x)
1250 double2Float (D# x) = F# (double2Float# x)
1251 float2Double (F# x) = D# (float2Double# x)
1252
1253 expDouble    (D# x) = D# (expDouble# x)
1254 logDouble    (D# x) = D# (logDouble# x)
1255 sqrtDouble   (D# x) = D# (sqrtDouble# x)
1256 sinDouble    (D# x) = D# (sinDouble# x)
1257 cosDouble    (D# x) = D# (cosDouble# x)
1258 tanDouble    (D# x) = D# (tanDouble# x)
1259 asinDouble   (D# x) = D# (asinDouble# x)
1260 acosDouble   (D# x) = D# (acosDouble# x)
1261 atanDouble   (D# x) = D# (atanDouble# x)
1262 sinhDouble   (D# x) = D# (sinhDouble# x)
1263 coshDouble   (D# x) = D# (coshDouble# x)
1264 tanhDouble   (D# x) = D# (tanhDouble# x)
1265
1266 powerDouble  (D# x) (D# y) = D# (x **## y)
1267 \end{code}