\begin{code}
#include "ieee-flpt.h"
module GHC.Float( module GHC.Float, Float(..), Double(..), Float#, Double#
, double2Int, int2Double, float2Int, int2Float )
where
import Data.Maybe
import Data.Bits
import GHC.Base
import GHC.List
import GHC.Enum
import GHC.Show
import GHC.Num
import GHC.Real
import GHC.Arr
import GHC.Float.RealFracMethods
import GHC.Float.ConversionUtils
import GHC.Integer.Logarithms ( integerLogBase# )
import GHC.Integer.Logarithms.Internals
infixr 8 **
\end{code}
%*********************************************************
%* *
\subsection{Standard numeric classes}
%* *
%*********************************************************
\begin{code}
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
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
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 :: a -> Bool
isInfinite :: a -> Bool
isDenormalized :: a -> Bool
isNegativeZero :: a -> Bool
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 0 x = x
scaleFloat k x
| isFix = x
| otherwise = encodeFloat m (n + clamp b k)
where (m,n) = decodeFloat x
(l,h) = floatRange x
d = floatDigits x
b = h l + 4*d
isFix = x == 0 || isNaN x || isInfinite x
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) ||
(x < 0 && isNegativeZero y) ||
(isNegativeZero x && isNegativeZero y)
= atan2 (y) x
| y == 0 && (x < 0 || isNegativeZero x)
= pi
| x==0 && y==0 = y
| otherwise = x + y
\end{code}
%*********************************************************
%* *
\subsection{Type @Float@}
%* *
%*********************************************************
\begin{code}
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
fromInteger i = F# (floatFromInteger i)
instance Real Float where
toRational (F# x#) =
case decodeFloat_Int# x# of
(# m#, e# #)
| e# >=# 0# ->
(smallInteger m# `shiftLInteger` e#) :% 1
| (int2Word# m# `and#` 1##) `eqWord#` 0## ->
case elimZerosInt# m# (negateInt# e#) of
(# n, d# #) -> n :% shiftLInteger 1 d#
| otherwise ->
smallInteger m# :% shiftLInteger 1 (negateInt# e#)
instance Fractional Float where
(/) x y = divideFloat x y
fromRational (n:%0)
| n == 0 = 0/0
| n < 0 = (1)/0
| otherwise = 1/0
fromRational (n:%d)
| n == 0 = encodeFloat 0 0
| n < 0 = (fromRat'' minEx mantDigs (n) d)
| otherwise = fromRat'' minEx mantDigs n d
where
minEx = FLT_MIN_EXP
mantDigs = FLT_MANT_DIG
recip x = 1.0 / x
instance RealFrac Float where
#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
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
asinh x = log (x + sqrt (1.0+x*x))
acosh x = log (x + (x+1.0) * sqrt ((x1.0)/(x+1.0)))
atanh x = 0.5 * log ((1.0+x) / (1.0x))
instance RealFloat Float where
floatRadix _ = FLT_RADIX
floatDigits _ = FLT_MANT_DIG
floatRange _ = (FLT_MIN_EXP, FLT_MAX_EXP)
decodeFloat (F# f#) = case decodeFloat_Int# f# of
(# i, e #) -> (smallInteger i, I# 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
significand x = case decodeFloat x of
(m,_) -> encodeFloat m (negate (floatDigits x))
scaleFloat 0 x = x
scaleFloat k x
| isFix = x
| otherwise = case decodeFloat x of
(m,n) -> encodeFloat m (n + clamp bf k)
where bf = FLT_MAX_EXP (FLT_MIN_EXP) + 4*FLT_MANT_DIG
isFix = x == 0 || isFloatFinite x == 0
isNaN x = 0 /= isFloatNaN x
isInfinite x = 0 /= isFloatInfinite x
isDenormalized x = 0 /= isFloatDenormalized x
isNegativeZero x = 0 /= isFloatNegativeZero x
isIEEE _ = True
instance Show Float where
showsPrec x = showSignedFloat showFloat x
showList = showList__ (showsPrec 0)
\end{code}
%*********************************************************
%* *
\subsection{Type @Double@}
%* *
%*********************************************************
\begin{code}
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
fromInteger i = D# (doubleFromInteger i)
instance Real Double where
toRational (D# x#) =
case decodeDoubleInteger x# of
(# m, e# #)
| e# >=# 0# ->
shiftLInteger m e# :% 1
| (integerToWord m `and#` 1##) `eqWord#` 0## ->
case elimZerosInteger m (negateInt# e#) of
(# n, d# #) -> n :% shiftLInteger 1 d#
| otherwise ->
m :% shiftLInteger 1 (negateInt# e#)
instance Fractional Double where
(/) x y = divideDouble x y
fromRational (n:%0)
| n == 0 = 0/0
| n < 0 = (1)/0
| otherwise = 1/0
fromRational (n:%d)
| n == 0 = encodeFloat 0 0
| n < 0 = (fromRat'' minEx mantDigs (n) d)
| otherwise = fromRat'' minEx mantDigs n d
where
minEx = DBL_MIN_EXP
mantDigs = DBL_MANT_DIG
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
asinh x = log (x + sqrt (1.0+x*x))
acosh x = log (x + (x+1.0) * sqrt ((x1.0)/(x+1.0)))
atanh x = 0.5 * log ((1.0+x) / (1.0x))
instance RealFrac Double where
properFraction x
= case (decodeFloat x) of { (m,n) ->
if n >= 0 then
(fromInteger m * 2 ^ n, 0.0)
else
case (quotRem m (2^(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
ceiling 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
floatDigits _ = DBL_MANT_DIG
floatRange _ = (DBL_MIN_EXP, DBL_MAX_EXP)
decodeFloat (D# x#)
= case decodeDoubleInteger x# of
(# i, j #) -> (i, I# j)
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
significand x = case decodeFloat x of
(m,_) -> encodeFloat m (negate (floatDigits x))
scaleFloat 0 x = x
scaleFloat k x
| isFix = x
| otherwise = case decodeFloat x of
(m,n) -> encodeFloat m (n + clamp bd k)
where bd = DBL_MAX_EXP (DBL_MIN_EXP) + 4*DBL_MANT_DIG
isFix = x == 0 || isDoubleFinite x == 0
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 = showSignedFloat showFloat x
showList = showList__ (showsPrec 0)
\end{code}
%*********************************************************
%* *
\subsection{@Enum@ instances}
%* *
%*********************************************************
The @Enum@ instances for Floats and Doubles are slightly unusual.
The @toEnum@ function truncates numbers to Int. The definitions
of @enumFrom@ and @enumFromThen@ allow floats to be used in arithmetic
series: [0,0.1 .. 1.0]. However, roundoff errors make these somewhat
dubious. This example may have either 10 or 11 elements, depending on
how 0.1 is represented.
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
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
toEnum = int2Float
fromEnum = fromInteger . truncate
enumFrom = numericEnumFrom
enumFromTo = numericEnumFromTo
enumFromThen = numericEnumFromThen
enumFromThenTo = numericEnumFromThenTo
instance Enum Double where
succ x = x + 1
pred x = x 1
toEnum = int2Double
fromEnum = fromInteger . truncate
enumFrom = numericEnumFrom
enumFromTo = numericEnumFromTo
enumFromThen = numericEnumFromThen
enumFromThenTo = numericEnumFromThenTo
\end{code}
%*********************************************************
%* *
\subsection{Printing floating point}
%* *
%*********************************************************
\begin{code}
showFloat :: (RealFloat a) => a -> ShowS
showFloat x = showString (formatRealFloat FFGeneric Nothing x)
data FFFormat = FFExponent | FFFixed | FFGeneric
formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
formatRealFloat fmt decs x
| 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
base = 10
doFmt format (is, e) =
let ds = map intToDigit is in
case format of
FFGeneric ->
doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
(is,e)
FFExponent ->
case decs of
Nothing ->
let show_e' = show (e1) in
case ds of
"0" -> "0.0e0"
[d] -> d : ".0e" ++ 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')
in
d:'.':ds' ++ 'e':show (e1+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 (n1) ('0':s) ""
f n s (r:rs) = f (n1) (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')
roundTo :: Int -> Int -> [Int] -> (Int,[Int])
roundTo base d is =
case f d True is of
x@(0,_) -> x
(1,xs) -> (1, 1:xs)
_ -> error "roundTo: bad Value"
where
b2 = base `quot` 2
f n _ [] = (0, replicate n 0)
f 0 e (x:xs) | x == b2 && e && all (== 0) xs = (0, [])
| otherwise = (if x >= b2 then 1 else 0, [])
f n _ (i:xs)
| i' == base = (1,0:ds)
| otherwise = (0,i':ds)
where
(c,ds) = f (n1) (even i) xs
i' = c + i
floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
floatToDigits _ 0 = ([0], 0)
floatToDigits base x =
let
(f0, e0) = decodeFloat x
(minExp0, _) = floatRange x
p = floatDigits x
b = floatRadix x
minExp = minExp0 p
(f, e) =
let n = minExp e0 in
if n > 0 then (f0 `quot` (expt b n), e0+n) else (f0, e0)
(r, s, mUp, mDn) =
if e >= 0 then
let be = expt b e in
if f == expt b (p1) then
(f*be*b*2, 2*b, be*b, be)
else
(f*be*2, 2, be, be)
else
if e > minExp && f == expt b (p1) then
(f*b*2, expt b (e+1)*2, b, 1)
else
(f*2, expt b (e)*2, 1, 1)
k :: Int
k =
let
k0 :: Int
k0 =
if b == 2 && base == 10 then
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) :: Float) +
fromIntegral e * log (fromInteger b)) /
log (fromInteger base))
fixup n =
if n >= 0 then
if r + mUp <= expt base n * s then n else fixup (n+1)
else
if expt base (n) * (r + mUp) <= s then n else fixup (n+1)
in
fixup k0
gen ds rn sN mUpN mDnN =
let
(dn, rn') = (rn * base) `quotRem` sN
mUpN' = mUpN * base
mDnN' = mDnN * base
in
case (rn' < mDnN', rn' + mUpN' > sN) of
(True, False) -> dn : ds
(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 =
if k >= 0 then
gen [] r (s * expt base k) mUp mDn
else
let bk = expt base (k) in
gen [] (r * bk) s (mUp * bk) (mDn * bk)
in
(map fromIntegral (reverse rds), k)
\end{code}
%*********************************************************
%* *
\subsection{Converting from a Rational to a RealFloat
%* *
%*********************************************************
[In response to a request for documentation of how fromRational works,
Joe Fasel writes:] A quite reasonable request! This code was added to
the Prelude just before the 1.2 release, when Lennart, working with an
early version of hbi, noticed that (read . show) was not the identity
for floating-point numbers. (There was a one-bit error about half the
time.) The original version of the conversion function was in fact
simply a floating-point divide, as you suggest above. The new version
is, I grant you, somewhat denser.
Unfortunately, Joe's code doesn't work! Here's an example:
main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
This program prints
0.0000000000000000
instead of
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))
\end{pseudocode}
Now, here's Lennart's code (which works)
\begin{code}
fromRat :: (RealFloat a) => Rational -> a
fromRat (n :% 0) | n > 0 = 1/0
| n < 0 = 1/0
| otherwise = 0/0
fromRat (n :% d) | n > 0 = fromRat' (n :% d)
| n < 0 = fromRat' ((n) :% d)
| otherwise = encodeFloat 0 0
fromRat' :: (RealFloat a) => Rational -> a
fromRat' x = r
where b = floatRadix r
p = floatDigits r
(minExp0, _) = floatRange r
minExp = minExp0 p
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
x0 = x / f
(x', p') = if x0 >= xMax then (x0 / toRational b, p0+1) else (x0, p0)
r = encodeFloat (round x') p'
minExpt, maxExpt :: Int
minExpt = 0
maxExpt = 1100
expt :: Integer -> Int -> Integer
expt base n =
if base == 2 && n >= minExpt && n <= maxExpt then
expts!n
else
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]]
integerLogBase :: Integer -> Integer -> Int
integerLogBase b i
| i < b = 0
| b == 2 = I# (integerLog2# i)
| otherwise = I# (integerLogBase# b i)
\end{code}
Unfortunately, the old conversion code was awfully slow due to
a) a slow integer logarithm
b) repeated calculation of gcd's
For the case of Rational's coming from a Float or Double via toRational,
we can exploit the fact that the denominator is a power of two, which for
these brings a huge speedup since we need only shift and add instead
of division.
The below is an adaption of fromRat' for the conversion to
Float or Double exploiting the known floatRadix and avoiding
divisions as much as possible.
\begin{code}
fromRat'' :: RealFloat a => Int -> Int -> Integer -> Integer -> a
fromRat'' minEx@(I# me#) mantDigs@(I# md#) n d =
case integerLog2IsPowerOf2# d of
(# ld#, pw# #)
| pw# ==# 0# ->
case integerLog2# n of
ln# | ln# >=# (ld# +# me# -# 1#) ->
if ln# <# md#
then encodeFloat n (I# (negateInt# ld#))
else let n' = n `shiftR` (I# (ln# +# 1# -# md#))
n'' = case roundingMode# n (ln# -# md#) of
0# -> n'
2# -> n' + 1
_ -> case fromInteger n' .&. (1 :: Int) of
0 -> n'
_ -> n' + 1
in encodeFloat n'' (I# (ln# -# ld# +# 1# -# md#))
| otherwise ->
case ld# +# (me# -# md#) of
ld'# | ld'# <=# 0# ->
encodeFloat n (I# ((me# -# md#) -# ld'#))
| ld'# <=# ln# ->
let n' = n `shiftR` (I# ld'#)
in case roundingMode# n (ld'# -# 1#) of
0# -> encodeFloat n' (minEx mantDigs)
1# -> if fromInteger n' .&. (1 :: Int) == 0
then encodeFloat n' (minExmantDigs)
else encodeFloat (n' + 1) (minExmantDigs)
_ -> encodeFloat (n' + 1) (minExmantDigs)
| ld'# ># (ln# +# 1#) -> encodeFloat 0 0
| otherwise ->
case integerLog2IsPowerOf2# n of
(# _, 0# #) -> encodeFloat 0 0
(# _, _ #) -> encodeFloat 1 (minEx mantDigs)
| otherwise ->
let ln = I# (integerLog2# n)
ld = I# ld#
p0 = max minEx (ln ld)
(n', d')
| p0 < mantDigs = (n `shiftL` (mantDigs p0), d)
| p0 == mantDigs = (n, d)
| otherwise = (n, d `shiftL` (p0 mantDigs))
scale p a b
| (b `shiftL` mantDigs) <= a = (p+1, a, b `shiftL` 1)
| otherwise = (p, a, b)
(p', n'', d'') = scale (p0mantDigs) n' d'
rdq = case n'' `quotRem` d'' of
(q,r) -> case compare (r `shiftL` 1) d'' of
LT -> q
EQ -> if fromInteger q .&. (1 :: Int) == 0
then q else q+1
GT -> q+1
in encodeFloat rdq p'
\end{code}
%*********************************************************
%* *
\subsection{Floating point numeric primops}
%* *
%*********************************************************
Definitions of the boxed PrimOps; these will be
used in the case of partial applications, etc.
\begin{code}
plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
plusFloat (F# x) (F# y) = F# (plusFloat# x y)
minusFloat (F# x) (F# y) = F# (minusFloat# x y)
timesFloat (F# x) (F# y) = F# (timesFloat# x y)
divideFloat (F# x) (F# y) = F# (divideFloat# x y)
negateFloat :: Float -> Float
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
expFloat, logFloat, sqrtFloat :: Float -> Float
sinFloat, cosFloat, tanFloat :: Float -> Float
asinFloat, acosFloat, atanFloat :: Float -> Float
sinhFloat, coshFloat, tanhFloat :: Float -> Float
expFloat (F# x) = F# (expFloat# x)
logFloat (F# x) = F# (logFloat# x)
sqrtFloat (F# x) = F# (sqrtFloat# x)
sinFloat (F# x) = F# (sinFloat# x)
cosFloat (F# x) = F# (cosFloat# x)
tanFloat (F# x) = F# (tanFloat# x)
asinFloat (F# x) = F# (asinFloat# x)
acosFloat (F# x) = F# (acosFloat# x)
atanFloat (F# x) = F# (atanFloat# x)
sinhFloat (F# x) = F# (sinhFloat# x)
coshFloat (F# x) = F# (coshFloat# x)
tanhFloat (F# x) = F# (tanhFloat# x)
powerFloat :: Float -> Float -> Float
powerFloat (F# x) (F# y) = F# (powerFloat# x y)
plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
plusDouble (D# x) (D# y) = D# (x +## y)
minusDouble (D# x) (D# y) = D# (x -## y)
timesDouble (D# x) (D# y) = D# (x *## y)
divideDouble (D# x) (D# y) = D# (x /## y)
negateDouble :: Double -> Double
negateDouble (D# x) = D# (negateDouble# x)
gtDouble, geDouble, eqDouble, neDouble, leDouble, ltDouble :: Double -> Double -> Bool
gtDouble (D# x) (D# y) = x >## y
geDouble (D# x) (D# y) = x >=## y
eqDouble (D# x) (D# y) = x ==## y
neDouble (D# x) (D# y) = x /=## y
ltDouble (D# x) (D# y) = x <## y
leDouble (D# x) (D# y) = x <=## y
double2Float :: Double -> Float
double2Float (D# x) = F# (double2Float# x)
float2Double :: Float -> Double
float2Double (F# x) = D# (float2Double# x)
expDouble, logDouble, sqrtDouble :: Double -> Double
sinDouble, cosDouble, tanDouble :: Double -> Double
asinDouble, acosDouble, atanDouble :: Double -> Double
sinhDouble, coshDouble, tanhDouble :: Double -> Double
expDouble (D# x) = D# (expDouble# x)
logDouble (D# x) = D# (logDouble# x)
sqrtDouble (D# x) = D# (sqrtDouble# x)
sinDouble (D# x) = D# (sinDouble# x)
cosDouble (D# x) = D# (cosDouble# x)
tanDouble (D# x) = D# (tanDouble# x)
asinDouble (D# x) = D# (asinDouble# x)
acosDouble (D# x) = D# (acosDouble# x)
atanDouble (D# x) = D# (atanDouble# x)
sinhDouble (D# x) = D# (sinhDouble# x)
coshDouble (D# x) = D# (coshDouble# x)
tanhDouble (D# x) = D# (tanhDouble# x)
powerDouble :: Double -> Double -> Double
powerDouble (D# x) (D# y) = D# (x **## y)
\end{code}
\begin{code}
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 "isFloatFinite" isFloatFinite :: Float -> Int
foreign import ccall unsafe "isDoubleNaN" isDoubleNaN :: Double -> Int
foreign import ccall unsafe "isDoubleInfinite" isDoubleInfinite :: Double -> Int
foreign import ccall unsafe "isDoubleDenormalized" isDoubleDenormalized :: Double -> Int
foreign import ccall unsafe "isDoubleNegativeZero" isDoubleNegativeZero :: Double -> Int
foreign import ccall unsafe "isDoubleFinite" isDoubleFinite :: Double -> Int
\end{code}
%*********************************************************
%* *
\subsection{Coercion rules}
%* *
%*********************************************************
\begin{code}
\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)
-> Int
-> a
-> 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}