ghc/lib/ghc/PrelNum.lhs

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