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 `quot` (expt b n), e0+n) else (f0, e0)
 599   (r, s, mUp, mDn) =
 600    if e >= 0 then
 601     let be = expt b e in
 602     if f == expt b (p-1) then
 603       (f*be*b*2, 2*b, be*b, be)     -- according to Burger and Dybvig
 604     else
 605       (f*be*2, 2, be, be)
 606    else
 607     if e > minExp && f == expt b (p-1) then
 608       (f*b*2, expt b (-e+1)*2, b, 1)
 609     else
 610       (f*2, expt 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 very slightly larger than 8651/28738
 618         -- (about 5.3558e-10), so if log x >= 0, the approximation
 619         -- k1 is too small, hence we add one and need one fixup step less.
 620         -- If log x < 0, the approximation errs rather on the high side.
 621         -- That is usually more than compensated for by ignoring the
 622         -- fractional part of logBase 2 x, but when x is a power of 1/2
 623         -- or slightly larger and the exponent is a multiple of the
 624         -- denominator of the rational approximation to logBase 10 2,
 625         -- k1 is larger than logBase 10 x. If k1 > 1 + logBase 10 x,
 626         -- we get a leading zero-digit we don't want.
 627         -- With the approximation 3/10, this happened for
 628         -- 0.5^1030, 0.5^1040, ..., 0.5^1070 and values close above.
 629         -- The approximation 8651/28738 guarantees k1 < 1 + logBase 10 x
 630         -- for IEEE-ish floating point types with exponent fields
 631         -- <= 17 bits and mantissae of several thousand bits, earlier
 632         -- convergents to logBase 10 2 would fail for long double.
 633         -- Using quot instead of div is a little faster and requires
 634         -- fewer fixup steps for negative lx.
 635         let lx = p - 1 + e0
 636             k1 = (lx * 8651) `quot` 28738
 637         in if lx >= 0 then k1 + 1 else k1
 638      else
 639         -- f :: Integer, log :: Float -> Float,
 640         --               ceiling :: Float -> Int
 641         ceiling ((log (fromInteger (f+1) :: Float) +
 642                  fromIntegral e * log (fromInteger b)) /
 643                    log (fromInteger base))
 644 --WAS:            fromInt e * log (fromInteger b))
 645
 646     fixup n =
 647       if n >= 0 then
 648         if r + mUp <= expt base n * s then n else fixup (n+1)
 649       else
 650         if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
 651    in
 652    fixup k0
 653
 654   gen ds rn sN mUpN mDnN =
 655    let
 656     (dn, rn') = (rn * base) `quotRem` sN
 657     mUpN' = mUpN * base
 658     mDnN' = mDnN * base
 659    in
 660    case (rn' < mDnN', rn' + mUpN' > sN) of
 661     (True,  False) -> dn : ds
 662     (False, True)  -> dn+1 : ds
 663     (True,  True)  -> if rn' * 2 < sN then dn : ds else dn+1 : ds
 664     (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
 665
 666   rds =
 667    if k >= 0 then
 668       gen [] r (s * expt base k) mUp mDn
 669    else
 670      let bk = expt base (-k) in
 671      gen [] (r * bk) s (mUp * bk) (mDn * bk)
 672  in
 673  (map fromIntegral (reverse rds), k)
 674
 675 \end{code}
 676
 677
 678 %*********************************************************
 679 %*                                                      *
 680 \subsection{Converting from a Rational to a RealFloat
 681 %*                                                      *
 682 %*********************************************************
 683
 684 [In response to a request for documentation of how fromRational works,
 685 Joe Fasel writes:] A quite reasonable request!  This code was added to
 686 the Prelude just before the 1.2 release, when Lennart, working with an
 687 early version of hbi, noticed that (read . show) was not the identity
 688 for floating-point numbers.  (There was a one-bit error about half the
 689 time.)  The original version of the conversion function was in fact
 690 simply a floating-point divide, as you suggest above. The new version
 691 is, I grant you, somewhat denser.
 692
 693 Unfortunately, Joe's code doesn't work!  Here's an example:
 694
 695 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
 696
 697 This program prints
 698         0.0000000000000000
 699 instead of
 700         1.8217369128763981e-300
 701
 702 Here's Joe's code:
 703
 704 \begin{pseudocode}
 705 fromRat :: (RealFloat a) => Rational -> a
 706 fromRat x = x'
 707         where x' = f e
 708
 709 --              If the exponent of the nearest floating-point number to x
 710 --              is e, then the significand is the integer nearest xb^(-e),
 711 --              where b is the floating-point radix.  We start with a good
 712 --              guess for e, and if it is correct, the exponent of the
 713 --              floating-point number we construct will again be e.  If
 714 --              not, one more iteration is needed.
 715
 716               f e   = if e' == e then y else f e'
 717                       where y      = encodeFloat (round (x * (1 % b)^^e)) e
 718                             (_,e') = decodeFloat y
 719               b     = floatRadix x'
 720
 721 --              We obtain a trial exponent by doing a floating-point
 722 --              division of x's numerator by its denominator.  The
 723 --              result of this division may not itself be the ultimate
 724 --              result, because of an accumulation of three rounding
 725 --              errors.
 726
 727               (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
 728                                         / fromInteger (denominator x))
 729 \end{pseudocode}
 730
 731 Now, here's Lennart's code (which works)
 732
 733 \begin{code}
 734 -- | Converts a 'Rational' value into any type in class 'RealFloat'.
 735 {-# SPECIALISE fromRat :: Rational -> Double,
 736                           Rational -> Float #-}
 737 fromRat :: (RealFloat a) => Rational -> a
 738
 739 -- Deal with special cases first, delegating the real work to fromRat'
 740 fromRat (n :% 0) | n > 0     =  1/0        -- +Infinity
 741                  | n < 0     = -1/0        -- -Infinity
 742                  | otherwise =  0/0        -- NaN
 743
 744 fromRat (n :% d) | n > 0     = fromRat' (n :% d)
 745                  | n < 0     = - fromRat' ((-n) :% d)
 746                  | otherwise = encodeFloat 0 0             -- Zero
 747
 748 -- Conversion process:
 749 -- Scale the rational number by the RealFloat base until
 750 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
 751 -- Then round the rational to an Integer and encode it with the exponent
 752 -- that we got from the scaling.
 753 -- To speed up the scaling process we compute the log2 of the number to get
 754 -- a first guess of the exponent.
 755
 756 fromRat' :: (RealFloat a) => Rational -> a
 757 -- Invariant: argument is strictly positive
 758 fromRat' x = r
 759   where b = floatRadix r
 760         p = floatDigits r
 761         (minExp0, _) = floatRange r
 762         minExp = minExp0 - p            -- the real minimum exponent
 763         xMin   = toRational (expt b (p-1))
 764         xMax   = toRational (expt b p)
 765         p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
 766         f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
 767         (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
 768         r = encodeFloat (round x') p'
 769
 770 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
 771 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
 772 scaleRat b minExp xMin xMax p x
 773  | p <= minExp = (x, p)
 774  | x >= xMax   = scaleRat b minExp xMin xMax (p+1) (x/b)
 775  | x < xMin    = scaleRat b minExp xMin xMax (p-1) (x*b)
 776  | otherwise   = (x, p)
 777
 778 -- Exponentiation with a cache for the most common numbers.
 779 minExpt, maxExpt :: Int
 780 minExpt = 0
 781 maxExpt = 1100
 782
 783 expt :: Integer -> Int -> Integer
 784 expt base n =
 785     if base == 2 && n >= minExpt && n <= maxExpt then
 786         expts!n
 787     else
 788         if base == 10 && n <= maxExpt10 then
 789             expts10!n
 790         else
 791             base^n
 792
 793 expts :: Array Int Integer
 794 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
 795
 796 maxExpt10 :: Int
 797 maxExpt10 = 324
 798
 799 expts10 :: Array Int Integer
 800 expts10 = array (minExpt,maxExpt10) [(n,10^n) | n <- [minExpt .. maxExpt10]]
 801
 802 -- Compute the (floor of the) log of i in base b.
 803 -- Simplest way would be just divide i by b until it's smaller then b, but that would
 804 -- be very slow!  We are just slightly more clever.
 805 integerLogBase :: Integer -> Integer -> Int
 806 integerLogBase b i
 807    | i < b     = 0
 808    | otherwise = doDiv (i `div` (b^l)) l
 809        where
 810         -- Try squaring the base first to cut down the number of divisions.
 811          l = 2 * integerLogBase (b*b) i
 812
 813          doDiv :: Integer -> Int -> Int
 814          doDiv x y
 815             | x < b     = y
 816             | otherwise = doDiv (x `div` b) (y+1)
 817
 818 \end{code}
 819
 820
 821 %*********************************************************
 822 %*                                                      *
 823 \subsection{Floating point numeric primops}
 824 %*                                                      *
 825 %*********************************************************
 826
 827 Definitions of the boxed PrimOps; these will be
 828 used in the case of partial applications, etc.
 829
 830 \begin{code}
 831 plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
 832 plusFloat   (F# x) (F# y) = F# (plusFloat# x y)
 833 minusFloat  (F# x) (F# y) = F# (minusFloat# x y)
 834 timesFloat  (F# x) (F# y) = F# (timesFloat# x y)
 835 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
 836
 837 negateFloat :: Float -> Float
 838 negateFloat (F# x)        = F# (negateFloat# x)
 839
 840 gtFloat, geFloat, eqFloat, neFloat, ltFloat, leFloat :: Float -> Float -> Bool
 841 gtFloat     (F# x) (F# y) = gtFloat# x y
 842 geFloat     (F# x) (F# y) = geFloat# x y
 843 eqFloat     (F# x) (F# y) = eqFloat# x y
 844 neFloat     (F# x) (F# y) = neFloat# x y
 845 ltFloat     (F# x) (F# y) = ltFloat# x y
 846 leFloat     (F# x) (F# y) = leFloat# x y
 847
 848 float2Int :: Float -> Int
 849 float2Int   (F# x) = I# (float2Int# x)
 850
 851 int2Float :: Int -> Float
 852 int2Float   (I# x) = F# (int2Float# x)
 853
 854 expFloat, logFloat, sqrtFloat :: Float -> Float
 855 sinFloat, cosFloat, tanFloat  :: Float -> Float
 856 asinFloat, acosFloat, atanFloat  :: Float -> Float
 857 sinhFloat, coshFloat, tanhFloat  :: Float -> Float
 858 expFloat    (F# x) = F# (expFloat# x)
 859 logFloat    (F# x) = F# (logFloat# x)
 860 sqrtFloat   (F# x) = F# (sqrtFloat# x)
 861 sinFloat    (F# x) = F# (sinFloat# x)
 862 cosFloat    (F# x) = F# (cosFloat# x)
 863 tanFloat    (F# x) = F# (tanFloat# x)
 864 asinFloat   (F# x) = F# (asinFloat# x)
 865 acosFloat   (F# x) = F# (acosFloat# x)
 866 atanFloat   (F# x) = F# (atanFloat# x)
 867 sinhFloat   (F# x) = F# (sinhFloat# x)
 868 coshFloat   (F# x) = F# (coshFloat# x)
 869 tanhFloat   (F# x) = F# (tanhFloat# x)
 870
 871 powerFloat :: Float -> Float -> Float
 872 powerFloat  (F# x) (F# y) = F# (powerFloat# x y)
 873
 874 -- definitions of the boxed PrimOps; these will be
 875 -- used in the case of partial applications, etc.
 876
 877 plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
 878 plusDouble   (D# x) (D# y) = D# (x +## y)
 879 minusDouble  (D# x) (D# y) = D# (x -## y)
 880 timesDouble  (D# x) (D# y) = D# (x *## y)
 881 divideDouble (D# x) (D# y) = D# (x /## y)
 882
 883 negateDouble :: Double -> Double
 884 negateDouble (D# x)        = D# (negateDouble# x)
 885
 886 gtDouble, geDouble, eqDouble, neDouble, leDouble, ltDouble :: Double -> Double -> Bool
 887 gtDouble    (D# x) (D# y) = x >## y
 888 geDouble    (D# x) (D# y) = x >=## y
 889 eqDouble    (D# x) (D# y) = x ==## y
 890 neDouble    (D# x) (D# y) = x /=## y
 891 ltDouble    (D# x) (D# y) = x <## y
 892 leDouble    (D# x) (D# y) = x <=## y
 893
 894 double2Int :: Double -> Int
 895 double2Int   (D# x) = I# (double2Int#   x)
 896
 897 int2Double :: Int -> Double
 898 int2Double   (I# x) = D# (int2Double#   x)
 899
 900 double2Float :: Double -> Float
 901 double2Float (D# x) = F# (double2Float# x)
 902
 903 float2Double :: Float -> Double
 904 float2Double (F# x) = D# (float2Double# x)
 905
 906 expDouble, logDouble, sqrtDouble :: Double -> Double
 907 sinDouble, cosDouble, tanDouble  :: Double -> Double
 908 asinDouble, acosDouble, atanDouble  :: Double -> Double
 909 sinhDouble, coshDouble, tanhDouble  :: Double -> Double
 910 expDouble    (D# x) = D# (expDouble# x)
 911 logDouble    (D# x) = D# (logDouble# x)
 912 sqrtDouble   (D# x) = D# (sqrtDouble# x)
 913 sinDouble    (D# x) = D# (sinDouble# x)
 914 cosDouble    (D# x) = D# (cosDouble# x)
 915 tanDouble    (D# x) = D# (tanDouble# x)
 916 asinDouble   (D# x) = D# (asinDouble# x)
 917 acosDouble   (D# x) = D# (acosDouble# x)
 918 atanDouble   (D# x) = D# (atanDouble# x)
 919 sinhDouble   (D# x) = D# (sinhDouble# x)
 920 coshDouble   (D# x) = D# (coshDouble# x)
 921 tanhDouble   (D# x) = D# (tanhDouble# x)
 922
 923 powerDouble :: Double -> Double -> Double
 924 powerDouble  (D# x) (D# y) = D# (x **## y)
 925 \end{code}
 926
 927 \begin{code}
 928 foreign import ccall unsafe "isFloatNaN" isFloatNaN :: Float -> Int
 929 foreign import ccall unsafe "isFloatInfinite" isFloatInfinite :: Float -> Int
 930 foreign import ccall unsafe "isFloatDenormalized" isFloatDenormalized :: Float -> Int
 931 foreign import ccall unsafe "isFloatNegativeZero" isFloatNegativeZero :: Float -> Int
 932
 933
 934 foreign import ccall unsafe "isDoubleNaN" isDoubleNaN :: Double -> Int
 935 foreign import ccall unsafe "isDoubleInfinite" isDoubleInfinite :: Double -> Int
 936 foreign import ccall unsafe "isDoubleDenormalized" isDoubleDenormalized :: Double -> Int
 937 foreign import ccall unsafe "isDoubleNegativeZero" isDoubleNegativeZero :: Double -> Int
 938 \end{code}
 939
 940 %*********************************************************
 941 %*                                                      *
 942 \subsection{Coercion rules}
 943 %*                                                      *
 944 %*********************************************************
 945
 946 \begin{code}
 947 {-# RULES
 948 "fromIntegral/Int->Float"   fromIntegral = int2Float
 949 "fromIntegral/Int->Double"  fromIntegral = int2Double
 950 "realToFrac/Float->Float"   realToFrac   = id :: Float -> Float
 951 "realToFrac/Float->Double"  realToFrac   = float2Double
 952 "realToFrac/Double->Float"  realToFrac   = double2Float
 953 "realToFrac/Double->Double" realToFrac   = id :: Double -> Double
 954 "realToFrac/Int->Double"    realToFrac   = int2Double   -- See Note [realToFrac int-to-float]
 955 "realToFrac/Int->Float"     realToFrac   = int2Float    --      ..ditto
 956     #-}
 957 \end{code}
 958
 959 Note [realToFrac int-to-float]
 960 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 961 Don found that the RULES for realToFrac/Int->Double and simliarly
 962 Float made a huge difference to some stream-fusion programs.  Here's
 963 an example
 964
 965       import Data.Array.Vector
 966
 967       n = 40000000
 968
 969       main = do
 970             let c = replicateU n (2::Double)
 971                 a = mapU realToFrac (enumFromToU 0 (n-1) ) :: UArr Double
 972             print (sumU (zipWithU (*) c a))
 973
 974 Without the RULE we get this loop body:
 975
 976       case $wtoRational sc_sY4 of ww_aM7 { (# ww1_aM9, ww2_aMa #) ->
 977       case $wfromRat ww1_aM9 ww2_aMa of tpl_X1P { D# ipv_sW3 ->
 978       Main.$s$wfold
 979         (+# sc_sY4 1)
 980         (+# wild_X1i 1)
 981         (+## sc2_sY6 (*## 2.0 ipv_sW3))
 982
 983 And with the rule:
 984
 985      Main.$s$wfold
 986         (+# sc_sXT 1)
 987         (+# wild_X1h 1)
 988         (+## sc2_sXV (*## 2.0 (int2Double# sc_sXT)))
 989
 990 The running time of the program goes from 120 seconds to 0.198 seconds
 991 with the native backend, and 0.143 seconds with the C backend.
 992
 993 A few more details in Trac #2251, and the patch message
 994 "Add RULES for realToFrac from Int".
 995
 996 %*********************************************************
 997 %*                                                      *
 998 \subsection{Utils}
 999 %*                                                      *
1000 %*********************************************************
1001
1002 \begin{code}
1003 showSignedFloat :: (RealFloat a)
1004   => (a -> ShowS)       -- ^ a function that can show unsigned values
1005   -> Int                -- ^ the precedence of the enclosing context
1006   -> a                  -- ^ the value to show
1007   -> ShowS
1008 showSignedFloat showPos p x
1009    | x < 0 || isNegativeZero x
1010        = showParen (p > 6) (showChar '-' . showPos (-x))
1011    | otherwise = showPos x
1012 \end{code}
1013
1014 We need to prevent over/underflow of the exponent in encodeFloat when
1015 called from scaleFloat, hence we clamp the scaling parameter.
1016 We must have a large enough range to cover the maximum difference of
1017 exponents returned by decodeFloat.
1018 \begin{code}
1019 clamp :: Int -> Int -> Int
1020 clamp bd k = max (-bd) (min bd k)
1021 \end{code}