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