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