GHC/Float.lhs

   1 \begin{code}
   2 {-# OPTIONS_GHC -XNoImplicitPrelude #-}
   3 -- We believe we could deorphan this module, by moving lots of things
   4 -- around, but we haven't got there yet:
   5 {-# OPTIONS_GHC -fno-warn-orphans #-}
   6 {-# OPTIONS_HADDOCK hide #-}
   7 -----------------------------------------------------------------------------
   8 -- |
   9 -- Module      :  GHC.Float
  10 -- Copyright   :  (c) The University of Glasgow 1994-2002
  11 -- License     :  see libraries/base/LICENSE
  12 --
  13 -- Maintainer  :  cvs-ghc@haskell.org
  14 -- Stability   :  internal
  15 -- Portability :  non-portable (GHC Extensions)
  16 --
  17 -- The types 'Float' and 'Double', and the classes 'Floating' and 'RealFloat'.
  18 --
  19 -----------------------------------------------------------------------------
  20
  21 #include "ieee-flpt.h"
  22
  23 -- #hide
  24 module GHC.Float( module GHC.Float, Float(..), Double(..), Float#, Double# )
  25     where
  26
  27 import Data.Maybe
  28
  29 import Data.Bits
  30 import GHC.Base
  31 import GHC.List
  32 import GHC.Enum
  33 import GHC.Show
  34 import GHC.Num
  35 import GHC.Real
  36 import GHC.Arr
  37
  38 infixr 8  **
  39 \end{code}
  40
  41 %*********************************************************
  42 %*                                                      *
  43 \subsection{Standard numeric classes}
  44 %*                                                      *
  45 %*********************************************************
  46
  47 \begin{code}
  48 -- | Trigonometric and hyperbolic functions and related functions.
  49 --
  50 -- Minimal complete definition:
  51 --      'pi', 'exp', 'log', 'sin', 'cos', 'sinh', 'cosh',
  52 --      'asin', 'acos', 'atan', 'asinh', 'acosh' and 'atanh'
  53 class  (Fractional a) => Floating a  where
  54     pi                  :: a
  55     exp, log, sqrt      :: a -> a
  56     (**), logBase       :: a -> a -> a
  57     sin, cos, tan       :: a -> a
  58     asin, acos, atan    :: a -> a
  59     sinh, cosh, tanh    :: a -> a
  60     asinh, acosh, atanh :: a -> a
  61
  62     {-# INLINE (**) #-}
  63     {-# INLINE logBase #-}
  64     {-# INLINE sqrt #-}
  65     {-# INLINE tan #-}
  66     {-# INLINE tanh #-}
  67     x ** y              =  exp (log x * y)
  68     logBase x y         =  log y / log x
  69     sqrt x              =  x ** 0.5
  70     tan  x              =  sin  x / cos  x
  71     tanh x              =  sinh x / cosh x
  72
  73 -- | Efficient, machine-independent access to the components of a
  74 -- floating-point number.
  75 --
  76 -- Minimal complete definition:
  77 --      all except 'exponent', 'significand', 'scaleFloat' and 'atan2'
  78 class  (RealFrac a, Floating a) => RealFloat a  where
  79     -- | a constant function, returning the radix of the representation
  80     -- (often @2@)
  81     floatRadix          :: a -> Integer
  82     -- | a constant function, returning the number of digits of
  83     -- 'floatRadix' in the significand
  84     floatDigits         :: a -> Int
  85     -- | a constant function, returning the lowest and highest values
  86     -- the exponent may assume
  87     floatRange          :: a -> (Int,Int)
  88     -- | The function 'decodeFloat' applied to a real floating-point
  89     -- number returns the significand expressed as an 'Integer' and an
  90     -- appropriately scaled exponent (an 'Int').  If @'decodeFloat' x@
  91     -- yields @(m,n)@, then @x@ is equal in value to @m*b^^n@, where @b@
  92     -- is the floating-point radix, and furthermore, either @m@ and @n@
  93     -- are both zero or else @b^(d-1) <= m < b^d@, where @d@ is the value
  94     -- of @'floatDigits' x@.  In particular, @'decodeFloat' 0 = (0,0)@.
  95     decodeFloat         :: a -> (Integer,Int)
  96     -- | 'encodeFloat' performs the inverse of 'decodeFloat'
  97     encodeFloat         :: Integer -> Int -> a
  98     -- | the second component of 'decodeFloat'.
  99     exponent            :: a -> Int
 100     -- | the first component of 'decodeFloat', scaled to lie in the open
 101     -- interval (@-1@,@1@)
 102     significand         :: a -> a
 103     -- | multiplies a floating-point number by an integer power of the radix
 104     scaleFloat          :: Int -> a -> a
 105     -- | 'True' if the argument is an IEEE \"not-a-number\" (NaN) value
 106     isNaN               :: a -> Bool
 107     -- | 'True' if the argument is an IEEE infinity or negative infinity
 108     isInfinite          :: a -> Bool
 109     -- | 'True' if the argument is too small to be represented in
 110     -- normalized format
 111     isDenormalized      :: a -> Bool
 112     -- | 'True' if the argument is an IEEE negative zero
 113     isNegativeZero      :: a -> Bool
 114     -- | 'True' if the argument is an IEEE floating point number
 115     isIEEE              :: a -> Bool
 116     -- | a version of arctangent taking two real floating-point arguments.
 117     -- For real floating @x@ and @y@, @'atan2' y x@ computes the angle
 118     -- (from the positive x-axis) of the vector from the origin to the
 119     -- point @(x,y)@.  @'atan2' y x@ returns a value in the range [@-pi@,
 120     -- @pi@].  It follows the Common Lisp semantics for the origin when
 121     -- signed zeroes are supported.  @'atan2' y 1@, with @y@ in a type
 122     -- that is 'RealFloat', should return the same value as @'atan' y@.
 123     -- A default definition of 'atan2' is provided, but implementors
 124     -- can provide a more accurate implementation.
 125     atan2               :: a -> a -> a
 126
 127
 128     exponent x          =  if m == 0 then 0 else n + floatDigits x
 129                            where (m,n) = decodeFloat x
 130
 131     significand x       =  encodeFloat m (negate (floatDigits x))
 132                            where (m,_) = decodeFloat x
 133
 134     scaleFloat k x      =  encodeFloat m (n + clamp b k)
 135                            where (m,n) = decodeFloat x
 136                                  (l,h) = floatRange x
 137                                  d     = floatDigits x
 138                                  b     = h - l + 4*d
 139                                  -- n+k may overflow, which would lead
 140                                  -- to wrong results, hence we clamp the
 141                                  -- scaling parameter.
 142                                  -- If n + k would be larger than h,
 143                                  -- n + clamp b k must be too, simliar
 144                                  -- for smaller than l - d.
 145                                  -- Add a little extra to keep clear
 146                                  -- from the boundary cases.
 147
 148     atan2 y x
 149       | x > 0            =  atan (y/x)
 150       | x == 0 && y > 0  =  pi/2
 151       | x <  0 && y > 0  =  pi + atan (y/x)
 152       |(x <= 0 && y < 0)            ||
 153        (x <  0 && isNegativeZero y) ||
 154        (isNegativeZero x && isNegativeZero y)
 155                          = -atan2 (-y) x
 156       | y == 0 && (x < 0 || isNegativeZero x)
 157                           =  pi    -- must be after the previous test on zero y
 158       | x==0 && y==0      =  y     -- must be after the other double zero tests
 159       | otherwise         =  x + y -- x or y is a NaN, return a NaN (via +)
 160 \end{code}
 161
 162
 163 %*********************************************************
 164 %*                                                      *
 165 \subsection{Type @Float@}
 166 %*                                                      *
 167 %*********************************************************
 168
 169 \begin{code}
 170 instance  Num Float  where
 171     (+)         x y     =  plusFloat x y
 172     (-)         x y     =  minusFloat x y
 173     negate      x       =  negateFloat x
 174     (*)         x y     =  timesFloat x y
 175     abs x | x >= 0.0    =  x
 176           | otherwise   =  negateFloat x
 177     signum x | x == 0.0  = 0
 178              | x > 0.0   = 1
 179              | otherwise = negate 1
 180
 181     {-# INLINE fromInteger #-}
 182     fromInteger i = F# (floatFromInteger i)
 183
 184 instance  Real Float  where
 185     toRational x        =  (m%1)*(b%1)^^n
 186                            where (m,n) = decodeFloat x
 187                                  b     = floatRadix  x
 188
 189 instance  Fractional Float  where
 190     (/) x y             =  divideFloat x y
 191     fromRational x      =  fromRat x
 192     recip x             =  1.0 / x
 193
 194 {-# RULES "truncate/Float->Int" truncate = float2Int #-}
 195 instance  RealFrac Float  where
 196
 197     {-# SPECIALIZE properFraction :: Float -> (Int, Float) #-}
 198     {-# SPECIALIZE round    :: Float -> Int #-}
 199
 200     {-# SPECIALIZE properFraction :: Float  -> (Integer, Float) #-}
 201     {-# SPECIALIZE round    :: Float -> Integer #-}
 202
 203         -- ceiling, floor, and truncate are all small
 204     {-# INLINE ceiling #-}
 205     {-# INLINE floor #-}
 206     {-# INLINE truncate #-}
 207
 208 -- We assume that FLT_RADIX is 2 so that we can use more efficient code
 209 #if FLT_RADIX != 2
 210 #error FLT_RADIX must be 2
 211 #endif
 212     properFraction (F# x#)
 213       = case decodeFloat_Int# x# of
 214         (# m#, n# #) ->
 215             let m = I# m#
 216                 n = I# n#
 217             in
 218             if n >= 0
 219             then (fromIntegral m * (2 ^ n), 0.0)
 220             else let i = if m >= 0 then                m `shiftR` negate n
 221                                    else negate (negate m `shiftR` negate n)
 222                      f = m - (i `shiftL` negate n)
 223                  in (fromIntegral i, encodeFloat (fromIntegral f) n)
 224
 225     truncate x  = case properFraction x of
 226                      (n,_) -> n
 227
 228     round x     = case properFraction x of
 229                      (n,r) -> let
 230                                 m         = if r < 0.0 then n - 1 else n + 1
 231                                 half_down = abs r - 0.5
 232                               in
 233                               case (compare half_down 0.0) of
 234                                 LT -> n
 235                                 EQ -> if even n then n else m
 236                                 GT -> m
 237
 238     ceiling x   = case properFraction x of
 239                     (n,r) -> if r > 0.0 then n + 1 else n
 240
 241     floor x     = case properFraction x of
 242                     (n,r) -> if r < 0.0 then n - 1 else n
 243
 244 instance  Floating Float  where
 245     pi                  =  3.141592653589793238
 246     exp x               =  expFloat x
 247     log x               =  logFloat x
 248     sqrt x              =  sqrtFloat x
 249     sin x               =  sinFloat x
 250     cos x               =  cosFloat x
 251     tan x               =  tanFloat x
 252     asin x              =  asinFloat x
 253     acos x              =  acosFloat x
 254     atan x              =  atanFloat x
 255     sinh x              =  sinhFloat x
 256     cosh x              =  coshFloat x
 257     tanh x              =  tanhFloat x
 258     (**) x y            =  powerFloat x y
 259     logBase x y         =  log y / log x
 260
 261     asinh x = log (x + sqrt (1.0+x*x))
 262     acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
 263     atanh x = 0.5 * log ((1.0+x) / (1.0-x))
 264
 265 instance  RealFloat Float  where
 266     floatRadix _        =  FLT_RADIX        -- from float.h
 267     floatDigits _       =  FLT_MANT_DIG     -- ditto
 268     floatRange _        =  (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
 269
 270     decodeFloat (F# f#) = case decodeFloat_Int# f# of
 271                           (# i, e #) -> (smallInteger i, I# e)
 272
 273     encodeFloat i (I# e) = F# (encodeFloatInteger i e)
 274
 275     exponent x          = case decodeFloat x of
 276                             (m,n) -> if m == 0 then 0 else n + floatDigits x
 277
 278     significand x       = case decodeFloat x of
 279                             (m,_) -> encodeFloat m (negate (floatDigits x))
 280
 281     scaleFloat k x      = case decodeFloat x of
 282                             (m,n) -> encodeFloat m (n + clamp bf k)
 283                         where bf = FLT_MAX_EXP - (FLT_MIN_EXP) + 4*FLT_MANT_DIG
 284
 285     isNaN x          = 0 /= isFloatNaN x
 286     isInfinite x     = 0 /= isFloatInfinite x
 287     isDenormalized x = 0 /= isFloatDenormalized x
 288     isNegativeZero x = 0 /= isFloatNegativeZero x
 289     isIEEE _         = True
 290
 291 instance  Show Float  where
 292     showsPrec   x = showSignedFloat showFloat x
 293     showList = showList__ (showsPrec 0)
 294 \end{code}
 295
 296 %*********************************************************
 297 %*                                                      *
 298 \subsection{Type @Double@}
 299 %*                                                      *
 300 %*********************************************************
 301
 302 \begin{code}
 303 instance  Num Double  where
 304     (+)         x y     =  plusDouble x y
 305     (-)         x y     =  minusDouble x y
 306     negate      x       =  negateDouble x
 307     (*)         x y     =  timesDouble x y
 308     abs x | x >= 0.0    =  x
 309           | otherwise   =  negateDouble x
 310     signum x | x == 0.0  = 0
 311              | x > 0.0   = 1
 312              | otherwise = negate 1
 313
 314     {-# INLINE fromInteger #-}
 315     fromInteger i = D# (doubleFromInteger i)
 316
 317
 318 instance  Real Double  where
 319     toRational x        =  (m%1)*(b%1)^^n
 320                            where (m,n) = decodeFloat x
 321                                  b     = floatRadix  x
 322
 323 instance  Fractional Double  where
 324     (/) x y             =  divideDouble x y
 325     fromRational x      =  fromRat x
 326     recip x             =  1.0 / x
 327
 328 instance  Floating Double  where
 329     pi                  =  3.141592653589793238
 330     exp x               =  expDouble x
 331     log x               =  logDouble x
 332     sqrt x              =  sqrtDouble x
 333     sin  x              =  sinDouble x
 334     cos  x              =  cosDouble x
 335     tan  x              =  tanDouble x
 336     asin x              =  asinDouble x
 337     acos x              =  acosDouble x
 338     atan x              =  atanDouble x
 339     sinh x              =  sinhDouble x
 340     cosh x              =  coshDouble x
 341     tanh x              =  tanhDouble x
 342     (**) x y            =  powerDouble x y
 343     logBase x y         =  log y / log x
 344
 345     asinh x = log (x + sqrt (1.0+x*x))
 346     acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
 347     atanh x = 0.5 * log ((1.0+x) / (1.0-x))
 348
 349 {-# RULES "truncate/Double->Int" truncate = double2Int #-}
 350 instance  RealFrac Double  where
 351
 352     {-# SPECIALIZE properFraction :: Double -> (Int, Double) #-}
 353     {-# SPECIALIZE round    :: Double -> Int #-}
 354
 355     {-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
 356     {-# SPECIALIZE round    :: Double -> Integer #-}
 357
 358         -- ceiling, floor, and truncate are all small
 359     {-# INLINE ceiling #-}
 360     {-# INLINE floor #-}
 361     {-# INLINE truncate #-}
 362
 363     properFraction x
 364       = case (decodeFloat x)      of { (m,n) ->
 365         let  b = floatRadix x     in
 366         if n >= 0 then
 367             (fromInteger m * fromInteger b ^ n, 0.0)
 368         else
 369             case (quotRem m (b^(negate n))) of { (w,r) ->
 370             (fromInteger w, encodeFloat r n)
 371             }
 372         }
 373
 374     truncate x  = case properFraction x of
 375                      (n,_) -> n
 376
 377     round x     = case properFraction x of
 378                      (n,r) -> let
 379                                 m         = if r < 0.0 then n - 1 else n + 1
 380                                 half_down = abs r - 0.5
 381                               in
 382                               case (compare half_down 0.0) of
 383                                 LT -> n
 384                                 EQ -> if even n then n else m
 385                                 GT -> m
 386
 387     ceiling x   = case properFraction x of
 388                     (n,r) -> if r > 0.0 then n + 1 else n
 389
 390     floor x     = case properFraction x of
 391                     (n,r) -> if r < 0.0 then n - 1 else n
 392
 393 instance  RealFloat Double  where
 394     floatRadix _        =  FLT_RADIX        -- from float.h
 395     floatDigits _       =  DBL_MANT_DIG     -- ditto
 396     floatRange _        =  (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
 397
 398     decodeFloat (D# x#)
 399       = case decodeDoubleInteger x#   of
 400           (# i, j #) -> (i, I# j)
 401
 402     encodeFloat i (I# j) = D# (encodeDoubleInteger i j)
 403
 404     exponent x          = case decodeFloat x of
 405                             (m,n) -> if m == 0 then 0 else n + floatDigits x
 406
 407     significand x       = case decodeFloat x of
 408                             (m,_) -> encodeFloat m (negate (floatDigits x))
 409
 410     scaleFloat k x      = case decodeFloat x of
 411                             (m,n) -> encodeFloat m (n + clamp bd k)
 412                         where bd = DBL_MAX_EXP - (DBL_MIN_EXP) + 4*DBL_MANT_DIG
 413
 414     isNaN x             = 0 /= isDoubleNaN x
 415     isInfinite x        = 0 /= isDoubleInfinite x
 416     isDenormalized x    = 0 /= isDoubleDenormalized x
 417     isNegativeZero x    = 0 /= isDoubleNegativeZero x
 418     isIEEE _            = True
 419
 420 instance  Show Double  where
 421     showsPrec   x = showSignedFloat showFloat x
 422     showList = showList__ (showsPrec 0)
 423 \end{code}
 424
 425 %*********************************************************
 426 %*                                                      *
 427 \subsection{@Enum@ instances}
 428 %*                                                      *
 429 %*********************************************************
 430
 431 The @Enum@ instances for Floats and Doubles are slightly unusual.
 432 The @toEnum@ function truncates numbers to Int.  The definitions
 433 of @enumFrom@ and @enumFromThen@ allow floats to be used in arithmetic
 434 series: [0,0.1 .. 1.0].  However, roundoff errors make these somewhat
 435 dubious.  This example may have either 10 or 11 elements, depending on
 436 how 0.1 is represented.
 437
 438 NOTE: The instances for Float and Double do not make use of the default
 439 methods for @enumFromTo@ and @enumFromThenTo@, as these rely on there being
 440 a `non-lossy' conversion to and from Ints. Instead we make use of the
 441 1.2 default methods (back in the days when Enum had Ord as a superclass)
 442 for these (@numericEnumFromTo@ and @numericEnumFromThenTo@ below.)
 443
 444 \begin{code}
 445 instance  Enum Float  where
 446     succ x         = x + 1
 447     pred x         = x - 1
 448     toEnum         = int2Float
 449     fromEnum       = fromInteger . truncate   -- may overflow
 450     enumFrom       = numericEnumFrom
 451     enumFromTo     = numericEnumFromTo
 452     enumFromThen   = numericEnumFromThen
 453     enumFromThenTo = numericEnumFromThenTo
 454
 455 instance  Enum Double  where
 456     succ x         = x + 1
 457     pred x         = x - 1
 458     toEnum         =  int2Double
 459     fromEnum       =  fromInteger . truncate   -- may overflow
 460     enumFrom       =  numericEnumFrom
 461     enumFromTo     =  numericEnumFromTo
 462     enumFromThen   =  numericEnumFromThen
 463     enumFromThenTo =  numericEnumFromThenTo
 464 \end{code}
 465
 466
 467 %*********************************************************
 468 %*                                                      *
 469 \subsection{Printing floating point}
 470 %*                                                      *
 471 %*********************************************************
 472
 473
 474 \begin{code}
 475 -- | Show a signed 'RealFloat' value to full precision
 476 -- using standard decimal notation for arguments whose absolute value lies
 477 -- between @0.1@ and @9,999,999@, and scientific notation otherwise.
 478 showFloat :: (RealFloat a) => a -> ShowS
 479 showFloat x  =  showString (formatRealFloat FFGeneric Nothing x)
 480
 481 -- These are the format types.  This type is not exported.
 482
 483 data FFFormat = FFExponent | FFFixed | FFGeneric
 484
 485 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
 486 formatRealFloat fmt decs x
 487    | isNaN x                   = "NaN"
 488    | isInfinite x              = if x < 0 then "-Infinity" else "Infinity"
 489    | x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
 490    | otherwise                 = doFmt fmt (floatToDigits (toInteger base) x)
 491  where
 492   base = 10
 493
 494   doFmt format (is, e) =
 495     let ds = map intToDigit is in
 496     case format of
 497      FFGeneric ->
 498       doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
 499             (is,e)
 500      FFExponent ->
 501       case decs of
 502        Nothing ->
 503         let show_e' = show (e-1) in
 504         case ds of
 505           "0"     -> "0.0e0"
 506           [d]     -> d : ".0e" ++ show_e'
 507           (d:ds') -> d : '.' : ds' ++ "e" ++ show_e'
 508           []      -> error "formatRealFloat/doFmt/FFExponent: []"
 509        Just dec ->
 510         let dec' = max dec 1 in
 511         case is of
 512          [0] -> '0' :'.' : take dec' (repeat '0') ++ "e0"
 513          _ ->
 514           let
 515            (ei,is') = roundTo base (dec'+1) is
 516            (d:ds') = map intToDigit (if ei > 0 then init is' else is')
 517           in
 518           d:'.':ds' ++ 'e':show (e-1+ei)
 519      FFFixed ->
 520       let
 521        mk0 ls = case ls of { "" -> "0" ; _ -> ls}
 522       in
 523       case decs of
 524        Nothing
 525           | e <= 0    -> "0." ++ replicate (-e) '0' ++ ds
 526           | otherwise ->
 527              let
 528                 f 0 s    rs  = mk0 (reverse s) ++ '.':mk0 rs
 529                 f n s    ""  = f (n-1) ('0':s) ""
 530                 f n s (r:rs) = f (n-1) (r:s) rs
 531              in
 532                 f e "" ds
 533        Just dec ->
 534         let dec' = max dec 0 in
 535         if e >= 0 then
 536          let
 537           (ei,is') = roundTo base (dec' + e) is
 538           (ls,rs)  = splitAt (e+ei) (map intToDigit is')
 539          in
 540          mk0 ls ++ (if null rs then "" else '.':rs)
 541         else
 542          let
 543           (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
 544           d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
 545          in
 546          d : (if null ds' then "" else '.':ds')
 547
 548
 549 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
 550 roundTo base d is =
 551   case f d is of
 552     x@(0,_) -> x
 553     (1,xs)  -> (1, 1:xs)
 554     _       -> error "roundTo: bad Value"
 555  where
 556   b2 = base `div` 2
 557
 558   f n []     = (0, replicate n 0)
 559   f 0 (x:_)  = (if x >= b2 then 1 else 0, [])
 560   f n (i:xs)
 561      | i' == base = (1,0:ds)
 562      | otherwise  = (0,i':ds)
 563       where
 564        (c,ds) = f (n-1) xs
 565        i'     = c + i
 566
 567 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
 568 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
 569 -- This version uses a much slower logarithm estimator. It should be improved.
 570
 571 -- | 'floatToDigits' takes a base and a non-negative 'RealFloat' number,
 572 -- and returns a list of digits and an exponent.
 573 -- In particular, if @x>=0@, and
 574 --
 575 -- > floatToDigits base x = ([d1,d2,...,dn], e)
 576 --
 577 -- then
 578 --
 579 --      (1) @n >= 1@
 580 --
 581 --      (2) @x = 0.d1d2...dn * (base**e)@
 582 --
 583 --      (3) @0 <= di <= base-1@
 584
 585 floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
 586 floatToDigits _ 0 = ([0], 0)
 587 floatToDigits base x =
 588  let
 589   (f0, e0) = decodeFloat x
 590   (minExp0, _) = floatRange x
 591   p = floatDigits x
 592   b = floatRadix x
 593   minExp = minExp0 - p -- the real minimum exponent
 594   -- Haskell requires that f be adjusted so denormalized numbers
 595   -- will have an impossibly low exponent.  Adjust for this.
 596   (f, e) =
 597    let n = minExp - e0 in
 598    if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
 599   (r, s, mUp, mDn) =
 600    if e >= 0 then
 601     let be = b^ e in
 602     if f == b^(p-1) then
 603       (f*be*b*2, 2*b, be*b, b)
 604     else
 605       (f*be*2, 2, be, be)
 606    else
 607     if e > minExp && f == b^(p-1) then
 608       (f*b*2, b^(-e+1)*2, b, 1)
 609     else
 610       (f*2, b^(-e)*2, 1, 1)
 611   k :: Int
 612   k =
 613    let
 614     k0 :: Int
 615     k0 =
 616      if b == 2 && base == 10 then
 617         -- logBase 10 2 is slightly bigger than 3/10 so
 618         -- the following will err on the low side.  Ignoring
 619         -- the fraction will make it err even more.
 620         -- Haskell promises that p-1 <= logBase b f < p.
 621         (p - 1 + e0) * 3 `div` 10
 622      else
 623         -- f :: Integer, log :: Float -> Float,
 624         --               ceiling :: Float -> Int
 625         ceiling ((log (fromInteger (f+1) :: Float) +
 626                  fromIntegral e * log (fromInteger b)) /
 627                    log (fromInteger base))
 628 --WAS:            fromInt e * log (fromInteger b))
 629
 630     fixup n =
 631       if n >= 0 then
 632         if r + mUp <= expt base n * s then n else fixup (n+1)
 633       else
 634         if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
 635    in
 636    fixup k0
 637
 638   gen ds rn sN mUpN mDnN =
 639    let
 640     (dn, rn') = (rn * base) `divMod` sN
 641     mUpN' = mUpN * base
 642     mDnN' = mDnN * base
 643    in
 644    case (rn' < mDnN', rn' + mUpN' > sN) of
 645     (True,  False) -> dn : ds
 646     (False, True)  -> dn+1 : ds
 647     (True,  True)  -> if rn' * 2 < sN then dn : ds else dn+1 : ds
 648     (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
 649
 650   rds =
 651    if k >= 0 then
 652       gen [] r (s * expt base k) mUp mDn
 653    else
 654      let bk = expt base (-k) in
 655      gen [] (r * bk) s (mUp * bk) (mDn * bk)
 656  in
 657  (map fromIntegral (reverse rds), k)
 658
 659 \end{code}
 660
 661
 662 %*********************************************************
 663 %*                                                      *
 664 \subsection{Converting from a Rational to a RealFloat
 665 %*                                                      *
 666 %*********************************************************
 667
 668 [In response to a request for documentation of how fromRational works,
 669 Joe Fasel writes:] A quite reasonable request!  This code was added to
 670 the Prelude just before the 1.2 release, when Lennart, working with an
 671 early version of hbi, noticed that (read . show) was not the identity
 672 for floating-point numbers.  (There was a one-bit error about half the
 673 time.)  The original version of the conversion function was in fact
 674 simply a floating-point divide, as you suggest above. The new version
 675 is, I grant you, somewhat denser.
 676
 677 Unfortunately, Joe's code doesn't work!  Here's an example:
 678
 679 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
 680
 681 This program prints
 682         0.0000000000000000
 683 instead of
 684         1.8217369128763981e-300
 685
 686 Here's Joe's code:
 687
 688 \begin{pseudocode}
 689 fromRat :: (RealFloat a) => Rational -> a
 690 fromRat x = x'
 691         where x' = f e
 692
 693 --              If the exponent of the nearest floating-point number to x
 694 --              is e, then the significand is the integer nearest xb^(-e),
 695 --              where b is the floating-point radix.  We start with a good
 696 --              guess for e, and if it is correct, the exponent of the
 697 --              floating-point number we construct will again be e.  If
 698 --              not, one more iteration is needed.
 699
 700               f e   = if e' == e then y else f e'
 701                       where y      = encodeFloat (round (x * (1 % b)^^e)) e
 702                             (_,e') = decodeFloat y
 703               b     = floatRadix x'
 704
 705 --              We obtain a trial exponent by doing a floating-point
 706 --              division of x's numerator by its denominator.  The
 707 --              result of this division may not itself be the ultimate
 708 --              result, because of an accumulation of three rounding
 709 --              errors.
 710
 711               (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
 712                                         / fromInteger (denominator x))
 713 \end{pseudocode}
 714
 715 Now, here's Lennart's code (which works)
 716
 717 \begin{code}
 718 -- | Converts a 'Rational' value into any type in class 'RealFloat'.
 719 {-# SPECIALISE fromRat :: Rational -> Double,
 720                           Rational -> Float #-}
 721 fromRat :: (RealFloat a) => Rational -> a
 722
 723 -- Deal with special cases first, delegating the real work to fromRat'
 724 fromRat (n :% 0) | n > 0     =  1/0        -- +Infinity
 725                  | n < 0     = -1/0        -- -Infinity
 726                  | otherwise =  0/0        -- NaN
 727
 728 fromRat (n :% d) | n > 0     = fromRat' (n :% d)
 729                  | n < 0     = - fromRat' ((-n) :% d)
 730                  | otherwise = encodeFloat 0 0             -- Zero
 731
 732 -- Conversion process:
 733 -- Scale the rational number by the RealFloat base until
 734 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
 735 -- Then round the rational to an Integer and encode it with the exponent
 736 -- that we got from the scaling.
 737 -- To speed up the scaling process we compute the log2 of the number to get
 738 -- a first guess of the exponent.
 739
 740 fromRat' :: (RealFloat a) => Rational -> a
 741 -- Invariant: argument is strictly positive
 742 fromRat' x = r
 743   where b = floatRadix r
 744         p = floatDigits r
 745         (minExp0, _) = floatRange r
 746         minExp = minExp0 - p            -- the real minimum exponent
 747         xMin   = toRational (expt b (p-1))
 748         xMax   = toRational (expt b p)
 749         p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
 750         f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
 751         (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
 752         r = encodeFloat (round x') p'
 753
 754 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
 755 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
 756 scaleRat b minExp xMin xMax p x
 757  | p <= minExp = (x, p)
 758  | x >= xMax   = scaleRat b minExp xMin xMax (p+1) (x/b)
 759  | x < xMin    = scaleRat b minExp xMin xMax (p-1) (x*b)
 760  | otherwise   = (x, p)
 761
 762 -- Exponentiation with a cache for the most common numbers.
 763 minExpt, maxExpt :: Int
 764 minExpt = 0
 765 maxExpt = 1100
 766
 767 expt :: Integer -> Int -> Integer
 768 expt base n =
 769     if base == 2 && n >= minExpt && n <= maxExpt then
 770         expts!n
 771     else
 772         base^n
 773
 774 expts :: Array Int Integer
 775 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
 776
 777 -- Compute the (floor of the) log of i in base b.
 778 -- Simplest way would be just divide i by b until it's smaller then b, but that would
 779 -- be very slow!  We are just slightly more clever.
 780 integerLogBase :: Integer -> Integer -> Int
 781 integerLogBase b i
 782    | i < b     = 0
 783    | otherwise = doDiv (i `div` (b^l)) l
 784        where
 785         -- Try squaring the base first to cut down the number of divisions.
 786          l = 2 * integerLogBase (b*b) i
 787
 788          doDiv :: Integer -> Int -> Int
 789          doDiv x y
 790             | x < b     = y
 791             | otherwise = doDiv (x `div` b) (y+1)
 792
 793 \end{code}
 794
 795
 796 %*********************************************************
 797 %*                                                      *
 798 \subsection{Floating point numeric primops}
 799 %*                                                      *
 800 %*********************************************************
 801
 802 Definitions of the boxed PrimOps; these will be
 803 used in the case of partial applications, etc.
 804
 805 \begin{code}
 806 plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
 807 plusFloat   (F# x) (F# y) = F# (plusFloat# x y)
 808 minusFloat  (F# x) (F# y) = F# (minusFloat# x y)
 809 timesFloat  (F# x) (F# y) = F# (timesFloat# x y)
 810 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
 811
 812 negateFloat :: Float -> Float
 813 negateFloat (F# x)        = F# (negateFloat# x)
 814
 815 gtFloat, geFloat, eqFloat, neFloat, ltFloat, leFloat :: Float -> Float -> Bool
 816 gtFloat     (F# x) (F# y) = gtFloat# x y
 817 geFloat     (F# x) (F# y) = geFloat# x y
 818 eqFloat     (F# x) (F# y) = eqFloat# x y
 819 neFloat     (F# x) (F# y) = neFloat# x y
 820 ltFloat     (F# x) (F# y) = ltFloat# x y
 821 leFloat     (F# x) (F# y) = leFloat# x y
 822
 823 float2Int :: Float -> Int
 824 float2Int   (F# x) = I# (float2Int# x)
 825
 826 int2Float :: Int -> Float
 827 int2Float   (I# x) = F# (int2Float# x)
 828
 829 expFloat, logFloat, sqrtFloat :: Float -> Float
 830 sinFloat, cosFloat, tanFloat  :: Float -> Float
 831 asinFloat, acosFloat, atanFloat  :: Float -> Float
 832 sinhFloat, coshFloat, tanhFloat  :: Float -> Float
 833 expFloat    (F# x) = F# (expFloat# x)
 834 logFloat    (F# x) = F# (logFloat# x)
 835 sqrtFloat   (F# x) = F# (sqrtFloat# x)
 836 sinFloat    (F# x) = F# (sinFloat# x)
 837 cosFloat    (F# x) = F# (cosFloat# x)
 838 tanFloat    (F# x) = F# (tanFloat# x)
 839 asinFloat   (F# x) = F# (asinFloat# x)
 840 acosFloat   (F# x) = F# (acosFloat# x)
 841 atanFloat   (F# x) = F# (atanFloat# x)
 842 sinhFloat   (F# x) = F# (sinhFloat# x)
 843 coshFloat   (F# x) = F# (coshFloat# x)
 844 tanhFloat   (F# x) = F# (tanhFloat# x)
 845
 846 powerFloat :: Float -> Float -> Float
 847 powerFloat  (F# x) (F# y) = F# (powerFloat# x y)
 848
 849 -- definitions of the boxed PrimOps; these will be
 850 -- used in the case of partial applications, etc.
 851
 852 plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
 853 plusDouble   (D# x) (D# y) = D# (x +## y)
 854 minusDouble  (D# x) (D# y) = D# (x -## y)
 855 timesDouble  (D# x) (D# y) = D# (x *## y)
 856 divideDouble (D# x) (D# y) = D# (x /## y)
 857
 858 negateDouble :: Double -> Double
 859 negateDouble (D# x)        = D# (negateDouble# x)
 860
 861 gtDouble, geDouble, eqDouble, neDouble, leDouble, ltDouble :: Double -> Double -> Bool
 862 gtDouble    (D# x) (D# y) = x >## y
 863 geDouble    (D# x) (D# y) = x >=## y
 864 eqDouble    (D# x) (D# y) = x ==## y
 865 neDouble    (D# x) (D# y) = x /=## y
 866 ltDouble    (D# x) (D# y) = x <## y
 867 leDouble    (D# x) (D# y) = x <=## y
 868
 869 double2Int :: Double -> Int
 870 double2Int   (D# x) = I# (double2Int#   x)
 871
 872 int2Double :: Int -> Double
 873 int2Double   (I# x) = D# (int2Double#   x)
 874
 875 double2Float :: Double -> Float
 876 double2Float (D# x) = F# (double2Float# x)
 877
 878 float2Double :: Float -> Double
 879 float2Double (F# x) = D# (float2Double# x)
 880
 881 expDouble, logDouble, sqrtDouble :: Double -> Double
 882 sinDouble, cosDouble, tanDouble  :: Double -> Double
 883 asinDouble, acosDouble, atanDouble  :: Double -> Double
 884 sinhDouble, coshDouble, tanhDouble  :: Double -> Double
 885 expDouble    (D# x) = D# (expDouble# x)
 886 logDouble    (D# x) = D# (logDouble# x)
 887 sqrtDouble   (D# x) = D# (sqrtDouble# x)
 888 sinDouble    (D# x) = D# (sinDouble# x)
 889 cosDouble    (D# x) = D# (cosDouble# x)
 890 tanDouble    (D# x) = D# (tanDouble# x)
 891 asinDouble   (D# x) = D# (asinDouble# x)
 892 acosDouble   (D# x) = D# (acosDouble# x)
 893 atanDouble   (D# x) = D# (atanDouble# x)
 894 sinhDouble   (D# x) = D# (sinhDouble# x)
 895 coshDouble   (D# x) = D# (coshDouble# x)
 896 tanhDouble   (D# x) = D# (tanhDouble# x)
 897
 898 powerDouble :: Double -> Double -> Double
 899 powerDouble  (D# x) (D# y) = D# (x **## y)
 900 \end{code}
 901
 902 \begin{code}
 903 foreign import ccall unsafe "isFloatNaN" isFloatNaN :: Float -> Int
 904 foreign import ccall unsafe "isFloatInfinite" isFloatInfinite :: Float -> Int
 905 foreign import ccall unsafe "isFloatDenormalized" isFloatDenormalized :: Float -> Int
 906 foreign import ccall unsafe "isFloatNegativeZero" isFloatNegativeZero :: Float -> Int
 907
 908
 909 foreign import ccall unsafe "isDoubleNaN" isDoubleNaN :: Double -> Int
 910 foreign import ccall unsafe "isDoubleInfinite" isDoubleInfinite :: Double -> Int
 911 foreign import ccall unsafe "isDoubleDenormalized" isDoubleDenormalized :: Double -> Int
 912 foreign import ccall unsafe "isDoubleNegativeZero" isDoubleNegativeZero :: Double -> Int
 913 \end{code}
 914
 915 %*********************************************************
 916 %*                                                      *
 917 \subsection{Coercion rules}
 918 %*                                                      *
 919 %*********************************************************
 920
 921 \begin{code}
 922 {-# RULES
 923 "fromIntegral/Int->Float"   fromIntegral = int2Float
 924 "fromIntegral/Int->Double"  fromIntegral = int2Double
 925 "realToFrac/Float->Float"   realToFrac   = id :: Float -> Float
 926 "realToFrac/Float->Double"  realToFrac   = float2Double
 927 "realToFrac/Double->Float"  realToFrac   = double2Float
 928 "realToFrac/Double->Double" realToFrac   = id :: Double -> Double
 929 "realToFrac/Int->Double"    realToFrac   = int2Double   -- See Note [realToFrac int-to-float]
 930 "realToFrac/Int->Float"     realToFrac   = int2Float    --      ..ditto
 931     #-}
 932 \end{code}
 933
 934 Note [realToFrac int-to-float]
 935 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 936 Don found that the RULES for realToFrac/Int->Double and simliarly
 937 Float made a huge difference to some stream-fusion programs.  Here's
 938 an example
 939
 940       import Data.Array.Vector
 941
 942       n = 40000000
 943
 944       main = do
 945             let c = replicateU n (2::Double)
 946                 a = mapU realToFrac (enumFromToU 0 (n-1) ) :: UArr Double
 947             print (sumU (zipWithU (*) c a))
 948
 949 Without the RULE we get this loop body:
 950
 951       case $wtoRational sc_sY4 of ww_aM7 { (# ww1_aM9, ww2_aMa #) ->
 952       case $wfromRat ww1_aM9 ww2_aMa of tpl_X1P { D# ipv_sW3 ->
 953       Main.$s$wfold
 954         (+# sc_sY4 1)
 955         (+# wild_X1i 1)
 956         (+## sc2_sY6 (*## 2.0 ipv_sW3))
 957
 958 And with the rule:
 959
 960      Main.$s$wfold
 961         (+# sc_sXT 1)
 962         (+# wild_X1h 1)
 963         (+## sc2_sXV (*## 2.0 (int2Double# sc_sXT)))
 964
 965 The running time of the program goes from 120 seconds to 0.198 seconds
 966 with the native backend, and 0.143 seconds with the C backend.
 967
 968 A few more details in Trac #2251, and the patch message
 969 "Add RULES for realToFrac from Int".
 970
 971 %*********************************************************
 972 %*                                                      *
 973 \subsection{Utils}
 974 %*                                                      *
 975 %*********************************************************
 976
 977 \begin{code}
 978 showSignedFloat :: (RealFloat a)
 979   => (a -> ShowS)       -- ^ a function that can show unsigned values
 980   -> Int                -- ^ the precedence of the enclosing context
 981   -> a                  -- ^ the value to show
 982   -> ShowS
 983 showSignedFloat showPos p x
 984    | x < 0 || isNegativeZero x
 985        = showParen (p > 6) (showChar '-' . showPos (-x))
 986    | otherwise = showPos x
 987 \end{code}
 988
 989 We need to prevent over/underflow of the exponent in encodeFloat when
 990 called from scaleFloat, hence we clamp the scaling parameter.
 991 We must have a large enough range to cover the maximum difference of
 992 exponents returned by decodeFloat.
 993 \begin{code}
 994 clamp :: Int -> Int -> Int
 995 clamp bd k = max (-bd) (min bd k)
 996 \end{code}