GHC/Float.lhs

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