GHC/Float.lhs

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