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