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