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