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