\begin{code}
-{-# OPTIONS -fno-implicit-prelude #-}
+{-# LANGUAGE CPP
+ , NoImplicitPrelude
+ , MagicHash
+ , UnboxedTuples
+ , ForeignFunctionInterface
+ #-}
+-- We believe we could deorphan this module, by moving lots of things
+-- around, but we haven't got there yet:
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+{-# OPTIONS_HADDOCK hide #-}
+
-----------------------------------------------------------------------------
-- |
-- Module : GHC.Float
-- Copyright : (c) The University of Glasgow 1994-2002
-- License : see libraries/base/LICENSE
---
+--
-- Maintainer : cvs-ghc@haskell.org
-- Stability : internal
-- Portability : non-portable (GHC Extensions)
#include "ieee-flpt.h"
-module GHC.Float( module GHC.Float, Float#, Double# ) where
+-- #hide
+module GHC.Float( module GHC.Float, Float(..), Double(..), Float#, Double# )
+ where
import Data.Maybe
+import Data.Bits
import GHC.Base
import GHC.List
import GHC.Enum
\end{code}
%*********************************************************
-%* *
+%* *
\subsection{Standard numeric classes}
-%* *
+%* *
%*********************************************************
\begin{code}
+-- | Trigonometric and hyperbolic functions and related functions.
+--
+-- Minimal complete definition:
+-- 'pi', 'exp', 'log', 'sin', 'cos', 'sinh', 'cosh',
+-- 'asin', 'acos', 'atan', 'asinh', 'acosh' and 'atanh'
class (Fractional a) => Floating a where
- pi :: a
- exp, log, sqrt :: a -> a
- (**), logBase :: a -> a -> a
- sin, cos, tan :: a -> a
- asin, acos, atan :: a -> a
- sinh, cosh, tanh :: a -> a
+ pi :: a
+ exp, log, sqrt :: a -> a
+ (**), logBase :: a -> a -> a
+ sin, cos, tan :: a -> a
+ asin, acos, atan :: a -> a
+ sinh, cosh, tanh :: a -> a
asinh, acosh, atanh :: a -> a
- x ** y = exp (log x * y)
- logBase x y = log y / log x
- sqrt x = x ** 0.5
- tan x = sin x / cos x
- tanh x = sinh x / cosh x
-
+ {-# INLINE (**) #-}
+ {-# INLINE logBase #-}
+ {-# INLINE sqrt #-}
+ {-# INLINE tan #-}
+ {-# INLINE tanh #-}
+ x ** y = exp (log x * y)
+ logBase x y = log y / log x
+ sqrt x = x ** 0.5
+ tan x = sin x / cos x
+ tanh x = sinh x / cosh x
+
+-- | Efficient, machine-independent access to the components of a
+-- floating-point number.
+--
+-- Minimal complete definition:
+-- all except 'exponent', 'significand', 'scaleFloat' and 'atan2'
class (RealFrac a, Floating a) => RealFloat a where
- floatRadix :: a -> Integer
- floatDigits :: a -> Int
- floatRange :: a -> (Int,Int)
- decodeFloat :: a -> (Integer,Int)
- encodeFloat :: Integer -> Int -> a
- exponent :: a -> Int
- significand :: a -> a
- scaleFloat :: Int -> a -> a
- isNaN, isInfinite, isDenormalized, isNegativeZero, isIEEE
- :: a -> Bool
- atan2 :: a -> a -> a
-
-
- exponent x = if m == 0 then 0 else n + floatDigits x
- where (m,n) = decodeFloat x
-
- significand x = encodeFloat m (negate (floatDigits x))
- where (m,_) = decodeFloat x
-
- scaleFloat k x = encodeFloat m (n+k)
- where (m,n) = decodeFloat x
-
+ -- | a constant function, returning the radix of the representation
+ -- (often @2@)
+ floatRadix :: a -> Integer
+ -- | a constant function, returning the number of digits of
+ -- 'floatRadix' in the significand
+ floatDigits :: a -> Int
+ -- | a constant function, returning the lowest and highest values
+ -- the exponent may assume
+ floatRange :: a -> (Int,Int)
+ -- | The function 'decodeFloat' applied to a real floating-point
+ -- number returns the significand expressed as an 'Integer' and an
+ -- appropriately scaled exponent (an 'Int'). If @'decodeFloat' x@
+ -- yields @(m,n)@, then @x@ is equal in value to @m*b^^n@, where @b@
+ -- is the floating-point radix, and furthermore, either @m@ and @n@
+ -- are both zero or else @b^(d-1) <= m < b^d@, where @d@ is the value
+ -- of @'floatDigits' x@. In particular, @'decodeFloat' 0 = (0,0)@.
+ decodeFloat :: a -> (Integer,Int)
+ -- | 'encodeFloat' performs the inverse of 'decodeFloat'
+ encodeFloat :: Integer -> Int -> a
+ -- | the second component of 'decodeFloat'.
+ exponent :: a -> Int
+ -- | the first component of 'decodeFloat', scaled to lie in the open
+ -- interval (@-1@,@1@)
+ significand :: a -> a
+ -- | multiplies a floating-point number by an integer power of the radix
+ scaleFloat :: Int -> a -> a
+ -- | 'True' if the argument is an IEEE \"not-a-number\" (NaN) value
+ isNaN :: a -> Bool
+ -- | 'True' if the argument is an IEEE infinity or negative infinity
+ isInfinite :: a -> Bool
+ -- | 'True' if the argument is too small to be represented in
+ -- normalized format
+ isDenormalized :: a -> Bool
+ -- | 'True' if the argument is an IEEE negative zero
+ isNegativeZero :: a -> Bool
+ -- | 'True' if the argument is an IEEE floating point number
+ isIEEE :: a -> Bool
+ -- | a version of arctangent taking two real floating-point arguments.
+ -- For real floating @x@ and @y@, @'atan2' y x@ computes the angle
+ -- (from the positive x-axis) of the vector from the origin to the
+ -- point @(x,y)@. @'atan2' y x@ returns a value in the range [@-pi@,
+ -- @pi@]. It follows the Common Lisp semantics for the origin when
+ -- signed zeroes are supported. @'atan2' y 1@, with @y@ in a type
+ -- that is 'RealFloat', should return the same value as @'atan' y@.
+ -- A default definition of 'atan2' is provided, but implementors
+ -- can provide a more accurate implementation.
+ atan2 :: a -> a -> a
+
+
+ exponent x = if m == 0 then 0 else n + floatDigits x
+ where (m,n) = decodeFloat x
+
+ significand x = encodeFloat m (negate (floatDigits x))
+ where (m,_) = decodeFloat x
+
+ scaleFloat k x = encodeFloat m (n + clamp b k)
+ where (m,n) = decodeFloat x
+ (l,h) = floatRange x
+ d = floatDigits x
+ b = h - l + 4*d
+ -- n+k may overflow, which would lead
+ -- to wrong results, hence we clamp the
+ -- scaling parameter.
+ -- If n + k would be larger than h,
+ -- n + clamp b k must be too, simliar
+ -- for smaller than l - d.
+ -- Add a little extra to keep clear
+ -- from the boundary cases.
+
atan2 y x
| x > 0 = atan (y/x)
| x == 0 && y > 0 = pi/2
- | x < 0 && y > 0 = pi + atan (y/x)
+ | x < 0 && y > 0 = pi + atan (y/x)
|(x <= 0 && y < 0) ||
(x < 0 && isNegativeZero y) ||
(isNegativeZero x && isNegativeZero y)
%*********************************************************
-%* *
-\subsection{Type @Integer@, @Float@, @Double@}
-%* *
-%*********************************************************
-
-\begin{code}
--- | Single-precision floating point numbers.
-data Float = F# Float#
-
--- | Double-precision floating point numbers.
-data Double = D# Double#
-
-instance CCallable Float
-instance CReturnable Float
-
-instance CCallable Double
-instance CReturnable Double
-\end{code}
-
-
-%*********************************************************
-%* *
+%* *
\subsection{Type @Float@}
-%* *
+%* *
%*********************************************************
\begin{code}
-instance Eq Float where
- (F# x) == (F# y) = x `eqFloat#` y
-
-instance Ord Float where
- (F# x) `compare` (F# y) | x `ltFloat#` y = LT
- | x `eqFloat#` y = EQ
- | otherwise = GT
-
- (F# x) < (F# y) = x `ltFloat#` y
- (F# x) <= (F# y) = x `leFloat#` y
- (F# x) >= (F# y) = x `geFloat#` y
- (F# x) > (F# y) = x `gtFloat#` y
-
instance Num Float where
- (+) x y = plusFloat x y
- (-) x y = minusFloat x y
- negate x = negateFloat x
- (*) x y = timesFloat x y
- abs x | x >= 0.0 = x
- | otherwise = negateFloat x
- signum x | x == 0.0 = 0
- | x > 0.0 = 1
- | otherwise = negate 1
+ (+) x y = plusFloat x y
+ (-) x y = minusFloat x y
+ negate x = negateFloat x
+ (*) x y = timesFloat x y
+ abs x | x >= 0.0 = x
+ | otherwise = negateFloat x
+ signum x | x == 0.0 = 0
+ | x > 0.0 = 1
+ | otherwise = negate 1
{-# INLINE fromInteger #-}
- fromInteger n = encodeFloat n 0
- -- It's important that encodeFloat inlines here, and that
- -- fromInteger in turn inlines,
- -- so that if fromInteger is applied to an (S# i) the right thing happens
+ fromInteger i = F# (floatFromInteger i)
instance Real Float where
- toRational x = (m%1)*(b%1)^^n
- where (m,n) = decodeFloat x
- b = floatRadix x
+ toRational x = (m%1)*(b%1)^^n
+ where (m,n) = decodeFloat x
+ b = floatRadix x
instance Fractional Float where
- (/) x y = divideFloat x y
- fromRational x = fromRat x
- recip x = 1.0 / x
+ (/) x y = divideFloat x y
+ fromRational x = fromRat x
+ recip x = 1.0 / x
{-# RULES "truncate/Float->Int" truncate = float2Int #-}
instance RealFrac Float where
{-# SPECIALIZE properFraction :: Float -> (Integer, Float) #-}
{-# SPECIALIZE round :: Float -> Integer #-}
- -- ceiling, floor, and truncate are all small
+ -- ceiling, floor, and truncate are all small
{-# INLINE ceiling #-}
{-# INLINE floor #-}
{-# INLINE truncate #-}
- properFraction x
- = case (decodeFloat x) of { (m,n) ->
- let b = floatRadix x in
- if n >= 0 then
- (fromInteger m * fromInteger b ^ n, 0.0)
- else
- case (quotRem m (b^(negate n))) of { (w,r) ->
- (fromInteger w, encodeFloat r n)
- }
- }
-
- truncate x = case properFraction x of
- (n,_) -> n
-
- round x = case properFraction x of
- (n,r) -> let
- m = if r < 0.0 then n - 1 else n + 1
- half_down = abs r - 0.5
- in
- case (compare half_down 0.0) of
- LT -> n
- EQ -> if even n then n else m
- GT -> m
+-- We assume that FLT_RADIX is 2 so that we can use more efficient code
+#if FLT_RADIX != 2
+#error FLT_RADIX must be 2
+#endif
+ properFraction (F# x#)
+ = case decodeFloat_Int# x# of
+ (# m#, n# #) ->
+ let m = I# m#
+ n = I# n#
+ in
+ if n >= 0
+ then (fromIntegral m * (2 ^ n), 0.0)
+ else let i = if m >= 0 then m `shiftR` negate n
+ else negate (negate m `shiftR` negate n)
+ f = m - (i `shiftL` negate n)
+ in (fromIntegral i, encodeFloat (fromIntegral f) n)
+
+ truncate x = case properFraction x of
+ (n,_) -> n
+
+ round x = case properFraction x of
+ (n,r) -> let
+ m = if r < 0.0 then n - 1 else n + 1
+ half_down = abs r - 0.5
+ in
+ case (compare half_down 0.0) of
+ LT -> n
+ EQ -> if even n then n else m
+ GT -> m
ceiling x = case properFraction x of
- (n,r) -> if r > 0.0 then n + 1 else n
+ (n,r) -> if r > 0.0 then n + 1 else n
- floor x = case properFraction x of
- (n,r) -> if r < 0.0 then n - 1 else n
+ floor x = case properFraction x of
+ (n,r) -> if r < 0.0 then n - 1 else n
instance Floating Float where
- pi = 3.141592653589793238
- exp x = expFloat x
- log x = logFloat x
- sqrt x = sqrtFloat x
- sin x = sinFloat x
- cos x = cosFloat x
- tan x = tanFloat x
- asin x = asinFloat x
- acos x = acosFloat x
- atan x = atanFloat x
- sinh x = sinhFloat x
- cosh x = coshFloat x
- tanh x = tanhFloat x
- (**) x y = powerFloat x y
- logBase x y = log y / log x
+ pi = 3.141592653589793238
+ exp x = expFloat x
+ log x = logFloat x
+ sqrt x = sqrtFloat x
+ sin x = sinFloat x
+ cos x = cosFloat x
+ tan x = tanFloat x
+ asin x = asinFloat x
+ acos x = acosFloat x
+ atan x = atanFloat x
+ sinh x = sinhFloat x
+ cosh x = coshFloat x
+ tanh x = tanhFloat x
+ (**) x y = powerFloat x y
+ logBase x y = log y / log x
asinh x = log (x + sqrt (1.0+x*x))
acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
- atanh x = log ((x+1.0) / sqrt (1.0-x*x))
+ atanh x = 0.5 * log ((1.0+x) / (1.0-x))
instance RealFloat Float where
- floatRadix _ = FLT_RADIX -- from float.h
- floatDigits _ = FLT_MANT_DIG -- ditto
- floatRange _ = (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
+ floatRadix _ = FLT_RADIX -- from float.h
+ floatDigits _ = FLT_MANT_DIG -- ditto
+ floatRange _ = (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
+
+ decodeFloat (F# f#) = case decodeFloat_Int# f# of
+ (# i, e #) -> (smallInteger i, I# e)
- decodeFloat (F# f#)
- = case decodeFloat# f# of
- (# exp#, s#, d# #) -> (J# s# d#, I# exp#)
+ encodeFloat i (I# e) = F# (encodeFloatInteger i e)
- encodeFloat (S# i) j = int_encodeFloat# i j
- encodeFloat (J# s# d#) e = encodeFloat# s# d# e
+ exponent x = case decodeFloat x of
+ (m,n) -> if m == 0 then 0 else n + floatDigits x
- exponent x = case decodeFloat x of
- (m,n) -> if m == 0 then 0 else n + floatDigits x
+ significand x = case decodeFloat x of
+ (m,_) -> encodeFloat m (negate (floatDigits x))
- significand x = case decodeFloat x of
- (m,_) -> encodeFloat m (negate (floatDigits x))
+ scaleFloat k x = case decodeFloat x of
+ (m,n) -> encodeFloat m (n + clamp bf k)
+ where bf = FLT_MAX_EXP - (FLT_MIN_EXP) + 4*FLT_MANT_DIG
- scaleFloat k x = case decodeFloat x of
- (m,n) -> encodeFloat m (n+k)
isNaN x = 0 /= isFloatNaN x
isInfinite x = 0 /= isFloatInfinite x
isDenormalized x = 0 /= isFloatDenormalized x
isIEEE _ = True
instance Show Float where
- showsPrec x = showSigned showFloat x
- showList = showList__ (showsPrec 0)
+ showsPrec x = showSignedFloat showFloat x
+ showList = showList__ (showsPrec 0)
\end{code}
%*********************************************************
-%* *
+%* *
\subsection{Type @Double@}
-%* *
+%* *
%*********************************************************
\begin{code}
-instance Eq Double where
- (D# x) == (D# y) = x ==## y
-
-instance Ord Double where
- (D# x) `compare` (D# y) | x <## y = LT
- | x ==## y = EQ
- | otherwise = GT
-
- (D# x) < (D# y) = x <## y
- (D# x) <= (D# y) = x <=## y
- (D# x) >= (D# y) = x >=## y
- (D# x) > (D# y) = x >## y
-
instance Num Double where
- (+) x y = plusDouble x y
- (-) x y = minusDouble x y
- negate x = negateDouble x
- (*) x y = timesDouble x y
- abs x | x >= 0.0 = x
- | otherwise = negateDouble x
- signum x | x == 0.0 = 0
- | x > 0.0 = 1
- | otherwise = negate 1
+ (+) x y = plusDouble x y
+ (-) x y = minusDouble x y
+ negate x = negateDouble x
+ (*) x y = timesDouble x y
+ abs x | x >= 0.0 = x
+ | otherwise = negateDouble x
+ signum x | x == 0.0 = 0
+ | x > 0.0 = 1
+ | otherwise = negate 1
{-# INLINE fromInteger #-}
- -- See comments with Num Float
- fromInteger (S# i#) = case (int2Double# i#) of { d# -> D# d# }
- fromInteger (J# s# d#) = encodeDouble# s# d# 0
+ fromInteger i = D# (doubleFromInteger i)
instance Real Double where
- toRational x = (m%1)*(b%1)^^n
- where (m,n) = decodeFloat x
- b = floatRadix x
+ toRational x = (m%1)*(b%1)^^n
+ where (m,n) = decodeFloat x
+ b = floatRadix x
instance Fractional Double where
- (/) x y = divideDouble x y
- fromRational x = fromRat x
- recip x = 1.0 / x
+ (/) x y = divideDouble x y
+ fromRational x = fromRat x
+ recip x = 1.0 / x
instance Floating Double where
- pi = 3.141592653589793238
- exp x = expDouble x
- log x = logDouble x
- sqrt x = sqrtDouble x
- sin x = sinDouble x
- cos x = cosDouble x
- tan x = tanDouble x
- asin x = asinDouble x
- acos x = acosDouble x
- atan x = atanDouble x
- sinh x = sinhDouble x
- cosh x = coshDouble x
- tanh x = tanhDouble x
- (**) x y = powerDouble x y
- logBase x y = log y / log x
+ pi = 3.141592653589793238
+ exp x = expDouble x
+ log x = logDouble x
+ sqrt x = sqrtDouble x
+ sin x = sinDouble x
+ cos x = cosDouble x
+ tan x = tanDouble x
+ asin x = asinDouble x
+ acos x = acosDouble x
+ atan x = atanDouble x
+ sinh x = sinhDouble x
+ cosh x = coshDouble x
+ tanh x = tanhDouble x
+ (**) x y = powerDouble x y
+ logBase x y = log y / log x
asinh x = log (x + sqrt (1.0+x*x))
acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
- atanh x = log ((x+1.0) / sqrt (1.0-x*x))
+ atanh x = 0.5 * log ((1.0+x) / (1.0-x))
{-# RULES "truncate/Double->Int" truncate = double2Int #-}
instance RealFrac Double where
{-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
{-# SPECIALIZE round :: Double -> Integer #-}
- -- ceiling, floor, and truncate are all small
+ -- ceiling, floor, and truncate are all small
{-# INLINE ceiling #-}
{-# INLINE floor #-}
{-# INLINE truncate #-}
properFraction x
= case (decodeFloat x) of { (m,n) ->
- let b = floatRadix x in
- if n >= 0 then
- (fromInteger m * fromInteger b ^ n, 0.0)
- else
- case (quotRem m (b^(negate n))) of { (w,r) ->
- (fromInteger w, encodeFloat r n)
- }
+ let b = floatRadix x in
+ if n >= 0 then
+ (fromInteger m * fromInteger b ^ n, 0.0)
+ else
+ case (quotRem m (b^(negate n))) of { (w,r) ->
+ (fromInteger w, encodeFloat r n)
+ }
}
- truncate x = case properFraction x of
- (n,_) -> n
+ truncate x = case properFraction x of
+ (n,_) -> n
- round x = case properFraction x of
- (n,r) -> let
- m = if r < 0.0 then n - 1 else n + 1
- half_down = abs r - 0.5
- in
- case (compare half_down 0.0) of
- LT -> n
- EQ -> if even n then n else m
- GT -> m
+ round x = case properFraction x of
+ (n,r) -> let
+ m = if r < 0.0 then n - 1 else n + 1
+ half_down = abs r - 0.5
+ in
+ case (compare half_down 0.0) of
+ LT -> n
+ EQ -> if even n then n else m
+ GT -> m
ceiling x = case properFraction x of
- (n,r) -> if r > 0.0 then n + 1 else n
+ (n,r) -> if r > 0.0 then n + 1 else n
- floor x = case properFraction x of
- (n,r) -> if r < 0.0 then n - 1 else n
+ floor x = case properFraction x of
+ (n,r) -> if r < 0.0 then n - 1 else n
instance RealFloat Double where
- floatRadix _ = FLT_RADIX -- from float.h
- floatDigits _ = DBL_MANT_DIG -- ditto
- floatRange _ = (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
+ floatRadix _ = FLT_RADIX -- from float.h
+ floatDigits _ = DBL_MANT_DIG -- ditto
+ floatRange _ = (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
decodeFloat (D# x#)
- = case decodeDouble# x# of
- (# exp#, s#, d# #) -> (J# s# d#, I# exp#)
+ = case decodeDoubleInteger x# of
+ (# i, j #) -> (i, I# j)
- encodeFloat (S# i) j = int_encodeDouble# i j
- encodeFloat (J# s# d#) e = encodeDouble# s# d# e
+ encodeFloat i (I# j) = D# (encodeDoubleInteger i j)
- exponent x = case decodeFloat x of
- (m,n) -> if m == 0 then 0 else n + floatDigits x
+ exponent x = case decodeFloat x of
+ (m,n) -> if m == 0 then 0 else n + floatDigits x
- significand x = case decodeFloat x of
- (m,_) -> encodeFloat m (negate (floatDigits x))
+ significand x = case decodeFloat x of
+ (m,_) -> encodeFloat m (negate (floatDigits x))
- scaleFloat k x = case decodeFloat x of
- (m,n) -> encodeFloat m (n+k)
+ scaleFloat k x = case decodeFloat x of
+ (m,n) -> encodeFloat m (n + clamp bd k)
+ where bd = DBL_MAX_EXP - (DBL_MIN_EXP) + 4*DBL_MANT_DIG
- isNaN x = 0 /= isDoubleNaN x
- isInfinite x = 0 /= isDoubleInfinite x
- isDenormalized x = 0 /= isDoubleDenormalized x
- isNegativeZero x = 0 /= isDoubleNegativeZero x
- isIEEE _ = True
+ isNaN x = 0 /= isDoubleNaN x
+ isInfinite x = 0 /= isDoubleInfinite x
+ isDenormalized x = 0 /= isDoubleDenormalized x
+ isNegativeZero x = 0 /= isDoubleNegativeZero x
+ isIEEE _ = True
instance Show Double where
- showsPrec x = showSigned showFloat x
- showList = showList__ (showsPrec 0)
+ showsPrec x = showSignedFloat showFloat x
+ showList = showList__ (showsPrec 0)
\end{code}
%*********************************************************
-%* *
+%* *
\subsection{@Enum@ instances}
-%* *
+%* *
%*********************************************************
The @Enum@ instances for Floats and Doubles are slightly unusual.
NOTE: The instances for Float and Double do not make use of the default
methods for @enumFromTo@ and @enumFromThenTo@, as these rely on there being
-a `non-lossy' conversion to and from Ints. Instead we make use of the
+a `non-lossy' conversion to and from Ints. Instead we make use of the
1.2 default methods (back in the days when Enum had Ord as a superclass)
for these (@numericEnumFromTo@ and @numericEnumFromThenTo@ below.)
\begin{code}
instance Enum Float where
- succ x = x + 1
- pred x = x - 1
+ succ x = x + 1
+ pred x = x - 1
toEnum = int2Float
fromEnum = fromInteger . truncate -- may overflow
- enumFrom = numericEnumFrom
+ enumFrom = numericEnumFrom
enumFromTo = numericEnumFromTo
enumFromThen = numericEnumFromThen
enumFromThenTo = numericEnumFromThenTo
instance Enum Double where
- succ x = x + 1
- pred x = x - 1
+ succ x = x + 1
+ pred x = x - 1
toEnum = int2Double
fromEnum = fromInteger . truncate -- may overflow
- enumFrom = numericEnumFrom
+ enumFrom = numericEnumFrom
enumFromTo = numericEnumFromTo
enumFromThen = numericEnumFromThen
enumFromThenTo = numericEnumFromThenTo
%*********************************************************
-%* *
+%* *
\subsection{Printing floating point}
-%* *
+%* *
%*********************************************************
\begin{code}
+-- | Show a signed 'RealFloat' value to full precision
+-- using standard decimal notation for arguments whose absolute value lies
+-- between @0.1@ and @9,999,999@, and scientific notation otherwise.
showFloat :: (RealFloat a) => a -> ShowS
showFloat x = showString (formatRealFloat FFGeneric Nothing x)
formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
formatRealFloat fmt decs x
- | isNaN x = "NaN"
+ | isNaN x = "NaN"
| isInfinite x = if x < 0 then "-Infinity" else "Infinity"
| x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
- | otherwise = doFmt fmt (floatToDigits (toInteger base) x)
- where
+ | otherwise = doFmt fmt (floatToDigits (toInteger base) x)
+ where
base = 10
doFmt format (is, e) =
case format of
FFGeneric ->
doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
- (is,e)
+ (is,e)
FFExponent ->
case decs of
Nothing ->
let show_e' = show (e-1) in
- case ds of
+ case ds of
"0" -> "0.0e0"
[d] -> d : ".0e" ++ show_e'
- (d:ds') -> d : '.' : ds' ++ "e" ++ show_e'
+ (d:ds') -> d : '.' : ds' ++ "e" ++ show_e'
+ [] -> error "formatRealFloat/doFmt/FFExponent: []"
Just dec ->
let dec' = max dec 1 in
case is of
[0] -> '0' :'.' : take dec' (repeat '0') ++ "e0"
_ ->
let
- (ei,is') = roundTo base (dec'+1) is
- (d:ds') = map intToDigit (if ei > 0 then init is' else is')
+ (ei,is') = roundTo base (dec'+1) is
+ (d:ds') = map intToDigit (if ei > 0 then init is' else is')
in
- d:'.':ds' ++ 'e':show (e-1+ei)
+ d:'.':ds' ++ 'e':show (e-1+ei)
FFFixed ->
let
mk0 ls = case ls of { "" -> "0" ; _ -> ls}
in
case decs of
Nothing
- | e <= 0 -> "0." ++ replicate (-e) '0' ++ ds
- | otherwise ->
- let
- f 0 s rs = mk0 (reverse s) ++ '.':mk0 rs
- f n s "" = f (n-1) ('0':s) ""
- f n s (r:rs) = f (n-1) (r:s) rs
- in
- f e "" ds
+ | e <= 0 -> "0." ++ replicate (-e) '0' ++ ds
+ | otherwise ->
+ let
+ f 0 s rs = mk0 (reverse s) ++ '.':mk0 rs
+ f n s "" = f (n-1) ('0':s) ""
+ f n s (r:rs) = f (n-1) (r:s) rs
+ in
+ f e "" ds
Just dec ->
let dec' = max dec 0 in
- if e >= 0 then
- let
- (ei,is') = roundTo base (dec' + e) is
- (ls,rs) = splitAt (e+ei) (map intToDigit is')
- in
- mk0 ls ++ (if null rs then "" else '.':rs)
- else
- let
- (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
- d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
- in
- d : (if null ds' then "" else '.':ds')
+ if e >= 0 then
+ let
+ (ei,is') = roundTo base (dec' + e) is
+ (ls,rs) = splitAt (e+ei) (map intToDigit is')
+ in
+ mk0 ls ++ (if null rs then "" else '.':rs)
+ else
+ let
+ (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
+ d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
+ in
+ d : (if null ds' then "" else '.':ds')
roundTo :: Int -> Int -> [Int] -> (Int,[Int])
case f d is of
x@(0,_) -> x
(1,xs) -> (1, 1:xs)
+ _ -> error "roundTo: bad Value"
where
b2 = base `div` 2
-- by R.G. Burger and R.K. Dybvig in PLDI 96.
-- This version uses a much slower logarithm estimator. It should be improved.
--- floatToDigits takes a base and a non-negative RealFloat number,
--- and returns a list of digits and an exponent.
--- In particular, if x>=0, and
--- floatToDigits base x = ([d1,d2,...,dn], e)
+-- | 'floatToDigits' takes a base and a non-negative 'RealFloat' number,
+-- and returns a list of digits and an exponent.
+-- In particular, if @x>=0@, and
+--
+-- > floatToDigits base x = ([d1,d2,...,dn], e)
+--
-- then
--- (a) n >= 1
--- (b) x = 0.d1d2...dn * (base**e)
--- (c) 0 <= di <= base-1
+--
+-- (1) @n >= 1@
+--
+-- (2) @x = 0.d1d2...dn * (base**e)@
+--
+-- (3) @0 <= di <= base-1@
floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
floatToDigits _ 0 = ([0], 0)
floatToDigits base x =
- let
+ let
(f0, e0) = decodeFloat x
(minExp0, _) = floatRange x
p = floatDigits x
minExp = minExp0 - p -- the real minimum exponent
-- Haskell requires that f be adjusted so denormalized numbers
-- will have an impossibly low exponent. Adjust for this.
- (f, e) =
+ (f, e) =
let n = minExp - e0 in
- if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
+ if n > 0 then (f0 `quot` (expt b n), e0+n) else (f0, e0)
(r, s, mUp, mDn) =
if e >= 0 then
- let be = b^ e in
- if f == b^(p-1) then
- (f*be*b*2, 2*b, be*b, b)
+ let be = expt b e in
+ if f == expt b (p-1) then
+ (f*be*b*2, 2*b, be*b, be) -- according to Burger and Dybvig
else
(f*be*2, 2, be, be)
else
- if e > minExp && f == b^(p-1) then
- (f*b*2, b^(-e+1)*2, b, 1)
+ if e > minExp && f == expt b (p-1) then
+ (f*b*2, expt b (-e+1)*2, b, 1)
else
- (f*2, b^(-e)*2, 1, 1)
+ (f*2, expt b (-e)*2, 1, 1)
+ k :: Int
k =
- let
+ let
+ k0 :: Int
k0 =
if b == 2 && base == 10 then
- -- logBase 10 2 is slightly bigger than 3/10 so
- -- the following will err on the low side. Ignoring
- -- the fraction will make it err even more.
- -- Haskell promises that p-1 <= logBase b f < p.
- (p - 1 + e0) * 3 `div` 10
+ -- logBase 10 2 is very slightly larger than 8651/28738
+ -- (about 5.3558e-10), so if log x >= 0, the approximation
+ -- k1 is too small, hence we add one and need one fixup step less.
+ -- If log x < 0, the approximation errs rather on the high side.
+ -- That is usually more than compensated for by ignoring the
+ -- fractional part of logBase 2 x, but when x is a power of 1/2
+ -- or slightly larger and the exponent is a multiple of the
+ -- denominator of the rational approximation to logBase 10 2,
+ -- k1 is larger than logBase 10 x. If k1 > 1 + logBase 10 x,
+ -- we get a leading zero-digit we don't want.
+ -- With the approximation 3/10, this happened for
+ -- 0.5^1030, 0.5^1040, ..., 0.5^1070 and values close above.
+ -- The approximation 8651/28738 guarantees k1 < 1 + logBase 10 x
+ -- for IEEE-ish floating point types with exponent fields
+ -- <= 17 bits and mantissae of several thousand bits, earlier
+ -- convergents to logBase 10 2 would fail for long double.
+ -- Using quot instead of div is a little faster and requires
+ -- fewer fixup steps for negative lx.
+ let lx = p - 1 + e0
+ k1 = (lx * 8651) `quot` 28738
+ in if lx >= 0 then k1 + 1 else k1
else
- ceiling ((log (fromInteger (f+1)) +
- fromInteger (int2Integer e) * log (fromInteger b)) /
- log (fromInteger base))
---WAS: fromInt e * log (fromInteger b))
+ -- f :: Integer, log :: Float -> Float,
+ -- ceiling :: Float -> Int
+ ceiling ((log (fromInteger (f+1) :: Float) +
+ fromIntegral e * log (fromInteger b)) /
+ log (fromInteger base))
+--WAS: fromInt e * log (fromInteger b))
fixup n =
if n >= 0 then
gen ds rn sN mUpN mDnN =
let
- (dn, rn') = (rn * base) `divMod` sN
+ (dn, rn') = (rn * base) `quotRem` sN
mUpN' = mUpN * base
mDnN' = mDnN * base
in
(False, True) -> dn+1 : ds
(True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
(False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
-
- rds =
+
+ rds =
if k >= 0 then
gen [] r (s * expt base k) mUp mDn
else
%*********************************************************
-%* *
+%* *
\subsection{Converting from a Rational to a RealFloat
-%* *
+%* *
%*********************************************************
[In response to a request for documentation of how fromRational works,
main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
This program prints
- 0.0000000000000000
+ 0.0000000000000000
instead of
- 1.8217369128763981e-300
+ 1.8217369128763981e-300
Here's Joe's code:
\begin{pseudocode}
fromRat :: (RealFloat a) => Rational -> a
fromRat x = x'
- where x' = f e
-
--- If the exponent of the nearest floating-point number to x
--- is e, then the significand is the integer nearest xb^(-e),
--- where b is the floating-point radix. We start with a good
--- guess for e, and if it is correct, the exponent of the
--- floating-point number we construct will again be e. If
--- not, one more iteration is needed.
-
- f e = if e' == e then y else f e'
- where y = encodeFloat (round (x * (1 % b)^^e)) e
- (_,e') = decodeFloat y
- b = floatRadix x'
-
--- We obtain a trial exponent by doing a floating-point
--- division of x's numerator by its denominator. The
--- result of this division may not itself be the ultimate
--- result, because of an accumulation of three rounding
--- errors.
-
- (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
- / fromInteger (denominator x))
+ where x' = f e
+
+-- If the exponent of the nearest floating-point number to x
+-- is e, then the significand is the integer nearest xb^(-e),
+-- where b is the floating-point radix. We start with a good
+-- guess for e, and if it is correct, the exponent of the
+-- floating-point number we construct will again be e. If
+-- not, one more iteration is needed.
+
+ f e = if e' == e then y else f e'
+ where y = encodeFloat (round (x * (1 % b)^^e)) e
+ (_,e') = decodeFloat y
+ b = floatRadix x'
+
+-- We obtain a trial exponent by doing a floating-point
+-- division of x's numerator by its denominator. The
+-- result of this division may not itself be the ultimate
+-- result, because of an accumulation of three rounding
+-- errors.
+
+ (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
+ / fromInteger (denominator x))
\end{pseudocode}
Now, here's Lennart's code (which works)
\begin{code}
+-- | Converts a 'Rational' value into any type in class 'RealFloat'.
{-# SPECIALISE fromRat :: Rational -> Double,
- Rational -> Float #-}
+ Rational -> Float #-}
fromRat :: (RealFloat a) => Rational -> a
-- Deal with special cases first, delegating the real work to fromRat'
-fromRat (n :% 0) | n > 0 = 1/0 -- +Infinity
- | n == 0 = 0/0 -- NaN
- | n < 0 = -1/0 -- -Infinity
+fromRat (n :% 0) | n > 0 = 1/0 -- +Infinity
+ | n < 0 = -1/0 -- -Infinity
+ | otherwise = 0/0 -- NaN
-fromRat (n :% d) | n > 0 = fromRat' (n :% d)
- | n == 0 = encodeFloat 0 0 -- Zero
- | n < 0 = - fromRat' ((-n) :% d)
+fromRat (n :% d) | n > 0 = fromRat' (n :% d)
+ | n < 0 = - fromRat' ((-n) :% d)
+ | otherwise = encodeFloat 0 0 -- Zero
-- Conversion process:
-- Scale the rational number by the RealFloat base until
fromRat' x = r
where b = floatRadix r
p = floatDigits r
- (minExp0, _) = floatRange r
- minExp = minExp0 - p -- the real minimum exponent
- xMin = toRational (expt b (p-1))
- xMax = toRational (expt b p)
- p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
- f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
- (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
- r = encodeFloat (round x') p'
+ (minExp0, _) = floatRange r
+ minExp = minExp0 - p -- the real minimum exponent
+ xMin = toRational (expt b (p-1))
+ xMax = toRational (expt b p)
+ p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
+ f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
+ (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
+ r = encodeFloat (round x') p'
-- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
-scaleRat b minExp xMin xMax p x
+scaleRat b minExp xMin xMax p x
| p <= minExp = (x, p)
| x >= xMax = scaleRat b minExp xMin xMax (p+1) (x/b)
| x < xMin = scaleRat b minExp xMin xMax (p-1) (x*b)
if base == 2 && n >= minExpt && n <= maxExpt then
expts!n
else
- base^n
+ if base == 10 && n <= maxExpt10 then
+ expts10!n
+ else
+ base^n
expts :: Array Int Integer
expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
+maxExpt10 :: Int
+maxExpt10 = 324
+
+expts10 :: Array Int Integer
+expts10 = array (minExpt,maxExpt10) [(n,10^n) | n <- [minExpt .. maxExpt10]]
+
-- Compute the (floor of the) log of i in base b.
-- Simplest way would be just divide i by b until it's smaller then b, but that would
-- be very slow! We are just slightly more clever.
| i < b = 0
| otherwise = doDiv (i `div` (b^l)) l
where
- -- Try squaring the base first to cut down the number of divisions.
+ -- Try squaring the base first to cut down the number of divisions.
l = 2 * integerLogBase (b*b) i
- doDiv :: Integer -> Int -> Int
- doDiv x y
- | x < b = y
- | otherwise = doDiv (x `div` b) (y+1)
+ doDiv :: Integer -> Int -> Int
+ doDiv x y
+ | x < b = y
+ | otherwise = doDiv (x `div` b) (y+1)
\end{code}
%*********************************************************
-%* *
+%* *
\subsection{Floating point numeric primops}
-%* *
+%* *
%*********************************************************
Definitions of the boxed PrimOps; these will be
negateFloat (F# x) = F# (negateFloat# x)
gtFloat, geFloat, eqFloat, neFloat, ltFloat, leFloat :: Float -> Float -> Bool
-gtFloat (F# x) (F# y) = gtFloat# x y
-geFloat (F# x) (F# y) = geFloat# x y
-eqFloat (F# x) (F# y) = eqFloat# x y
-neFloat (F# x) (F# y) = neFloat# x y
-ltFloat (F# x) (F# y) = ltFloat# x y
-leFloat (F# x) (F# y) = leFloat# x y
+gtFloat (F# x) (F# y) = gtFloat# x y
+geFloat (F# x) (F# y) = geFloat# x y
+eqFloat (F# x) (F# y) = eqFloat# x y
+neFloat (F# x) (F# y) = neFloat# x y
+ltFloat (F# x) (F# y) = ltFloat# x y
+leFloat (F# x) (F# y) = leFloat# x y
float2Int :: Float -> Int
float2Int (F# x) = I# (float2Int# x)
\end{code}
\begin{code}
-foreign import ccall unsafe "__encodeFloat"
- encodeFloat# :: Int# -> ByteArray# -> Int -> Float
-foreign import ccall unsafe "__int_encodeFloat"
- int_encodeFloat# :: Int# -> Int -> Float
-
-
foreign import ccall unsafe "isFloatNaN" isFloatNaN :: Float -> Int
foreign import ccall unsafe "isFloatInfinite" isFloatInfinite :: Float -> Int
foreign import ccall unsafe "isFloatDenormalized" isFloatDenormalized :: Float -> Int
foreign import ccall unsafe "isFloatNegativeZero" isFloatNegativeZero :: Float -> Int
-foreign import ccall unsafe "__encodeDouble"
- encodeDouble# :: Int# -> ByteArray# -> Int -> Double
-foreign import ccall unsafe "__int_encodeDouble"
- int_encodeDouble# :: Int# -> Int -> Double
-
foreign import ccall unsafe "isDoubleNaN" isDoubleNaN :: Double -> Int
foreign import ccall unsafe "isDoubleInfinite" isDoubleInfinite :: Double -> Int
foreign import ccall unsafe "isDoubleDenormalized" isDoubleDenormalized :: Double -> Int
\end{code}
%*********************************************************
-%* *
+%* *
\subsection{Coercion rules}
-%* *
+%* *
%*********************************************************
\begin{code}
"realToFrac/Float->Double" realToFrac = float2Double
"realToFrac/Double->Float" realToFrac = double2Float
"realToFrac/Double->Double" realToFrac = id :: Double -> Double
+"realToFrac/Int->Double" realToFrac = int2Double -- See Note [realToFrac int-to-float]
+"realToFrac/Int->Float" realToFrac = int2Float -- ..ditto
#-}
\end{code}
+
+Note [realToFrac int-to-float]
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Don found that the RULES for realToFrac/Int->Double and simliarly
+Float made a huge difference to some stream-fusion programs. Here's
+an example
+
+ import Data.Array.Vector
+
+ n = 40000000
+
+ main = do
+ let c = replicateU n (2::Double)
+ a = mapU realToFrac (enumFromToU 0 (n-1) ) :: UArr Double
+ print (sumU (zipWithU (*) c a))
+
+Without the RULE we get this loop body:
+
+ case $wtoRational sc_sY4 of ww_aM7 { (# ww1_aM9, ww2_aMa #) ->
+ case $wfromRat ww1_aM9 ww2_aMa of tpl_X1P { D# ipv_sW3 ->
+ Main.$s$wfold
+ (+# sc_sY4 1)
+ (+# wild_X1i 1)
+ (+## sc2_sY6 (*## 2.0 ipv_sW3))
+
+And with the rule:
+
+ Main.$s$wfold
+ (+# sc_sXT 1)
+ (+# wild_X1h 1)
+ (+## sc2_sXV (*## 2.0 (int2Double# sc_sXT)))
+
+The running time of the program goes from 120 seconds to 0.198 seconds
+with the native backend, and 0.143 seconds with the C backend.
+
+A few more details in Trac #2251, and the patch message
+"Add RULES for realToFrac from Int".
+
+%*********************************************************
+%* *
+\subsection{Utils}
+%* *
+%*********************************************************
+
+\begin{code}
+showSignedFloat :: (RealFloat a)
+ => (a -> ShowS) -- ^ a function that can show unsigned values
+ -> Int -- ^ the precedence of the enclosing context
+ -> a -- ^ the value to show
+ -> ShowS
+showSignedFloat showPos p x
+ | x < 0 || isNegativeZero x
+ = showParen (p > 6) (showChar '-' . showPos (-x))
+ | otherwise = showPos x
+\end{code}
+
+We need to prevent over/underflow of the exponent in encodeFloat when
+called from scaleFloat, hence we clamp the scaling parameter.
+We must have a large enough range to cover the maximum difference of
+exponents returned by decodeFloat.
+\begin{code}
+clamp :: Int -> Int -> Int
+clamp bd k = max (-bd) (min bd k)
+\end{code}
projects / ghc-base.git / blobdiff