\begin{code}
-{-# OPTIONS_GHC -fno-implicit-prelude #-}
+{-# OPTIONS_GHC -XNoImplicitPrelude #-}
+-- 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
#include "ieee-flpt.h"
-- #hide
-module GHC.Float( module GHC.Float, Float#, Double# ) where
+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'
+-- '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'
+-- all except 'exponent', 'significand', 'scaleFloat' and 'atan2'
class (RealFrac a, Floating a) => RealFloat a where
-- | a constant function, returning the radix of the representation
-- (often @2@)
- floatRadix :: a -> Integer
+ floatRadix :: a -> Integer
-- | a constant function, returning the number of digits of
-- 'floatRadix' in the significand
- floatDigits :: a -> Int
+ floatDigits :: a -> Int
-- | a constant function, returning the lowest and highest values
-- the exponent may assume
- floatRange :: a -> (Int,Int)
+ 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@
-- 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)
+ decodeFloat :: a -> (Integer,Int)
-- | 'encodeFloat' performs the inverse of 'decodeFloat'
- encodeFloat :: Integer -> Int -> a
+ encodeFloat :: Integer -> Int -> a
-- | the second component of 'decodeFloat'.
- exponent :: a -> Int
+ exponent :: a -> Int
-- | the first component of 'decodeFloat', scaled to lie in the open
-- interval (@-1@,@1@)
- significand :: a -> a
+ significand :: a -> a
-- | multiplies a floating-point number by an integer power of the radix
- scaleFloat :: Int -> a -> a
+ scaleFloat :: Int -> a -> a
-- | 'True' if the argument is an IEEE \"not-a-number\" (NaN) value
- isNaN :: a -> Bool
+ isNaN :: a -> Bool
-- | 'True' if the argument is an IEEE infinity or negative infinity
- isInfinite :: a -> Bool
+ isInfinite :: a -> Bool
-- | 'True' if the argument is too small to be represented in
-- normalized format
- isDenormalized :: a -> Bool
+ isDenormalized :: a -> Bool
-- | 'True' if the argument is an IEEE negative zero
- isNegativeZero :: a -> Bool
+ isNegativeZero :: a -> Bool
-- | 'True' if the argument is an IEEE floating point number
- isIEEE :: a -> Bool
+ 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
-- 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
+ atan2 :: a -> a -> a
- exponent x = if m == 0 then 0 else n + floatDigits x
- where (m,n) = decodeFloat x
+ 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
+ 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.
- scaleFloat k x = encodeFloat m (n+k)
- where (m,n) = decodeFloat x
-
atan2 y x
| x > 0 = atan (y/x)
| x == 0 && y > 0 = pi/2
%*********************************************************
-%* *
-\subsection{Type @Integer@, @Float@, @Double@}
-%* *
-%*********************************************************
-
-\begin{code}
--- | Single-precision floating point numbers.
--- It is desirable that this type be at least equal in range and precision
--- to the IEEE single-precision type.
-data Float = F# Float#
-
--- | Double-precision floating point numbers.
--- It is desirable that this type be at least equal in range and precision
--- to the IEEE double-precision type.
-data Double = D# 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# f# of
- (# exp#, s#, d# #) -> (J# s# d#, I# exp#)
+ decodeFloat (F# f#) = case decodeFloat_Int# f# of
+ (# i, e #) -> (smallInteger i, I# e)
- encodeFloat (S# i) j = int_encodeFloat# i j
- encodeFloat (J# s# d#) e = encodeFloat# s# d# e
+ encodeFloat i (I# e) = F# (encodeFloatInteger i 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
+ 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
+ showsPrec x = showSignedFloat showFloat x
showList = showList__ (showsPrec 0)
\end{code}
%*********************************************************
-%* *
+%* *
\subsection{@Enum@ instances}
-%* *
+%* *
%*********************************************************
The @Enum@ instances for Floats and Doubles are slightly unusual.
\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}
-%* *
+%* *
%*********************************************************
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)
+ | otherwise = doFmt fmt (floatToDigits (toInteger base) x)
where
base = 10
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
--
-- then
--
--- (1) @n >= 1@
+-- (1) @n >= 1@
--
--- (2) @x = 0.d1d2...dn * (base**e)@
+-- (2) @x = 0.d1d2...dn * (base**e)@
--
--- (3) @0 <= di <= base-1@
+-- (3) @0 <= di <= base-1@
floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
floatToDigits _ 0 = ([0], 0)
(f*b*2, b^(-e+1)*2, b, 1)
else
(f*2, b^(-e)*2, 1, 1)
+ k :: Int
k =
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
+ -- 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
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
%*********************************************************
-%* *
+%* *
\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)
| 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