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