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