\begin{code}
-{-# OPTIONS_GHC -XNoImplicitPrelude #-}
+{-# 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"
-- #hide
-module GHC.Float( module GHC.Float, Float(..), Double(..), Float#, Double# )
+module GHC.Float( module GHC.Float, Float(..), Double(..), Float#, Double#
+ , double2Int, int2Double, float2Int, int2Float )
where
import Data.Maybe
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}
significand x = encodeFloat m (negate (floatDigits x))
where (m,_) = decodeFloat x
- scaleFloat k x = encodeFloat m (n+k)
+ 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)
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 (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 x = fromRat x
+ 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
-{-# RULES "truncate/Float->Int" truncate = float2Int #-}
+-- RULES for Integer and Int
+{-# RULES
+"properFraction/Float->Integer" properFraction = properFractionFloatInteger
+"truncate/Float->Integer" truncate = truncateFloatInteger
+"floor/Float->Integer" floor = floorFloatInteger
+"ceiling/Float->Integer" ceiling = ceilingFloatInteger
+"round/Float->Integer" round = roundFloatInteger
+"properFraction/Float->Int" properFraction = properFractionFloatInt
+"truncate/Float->Int" truncate = float2Int
+"floor/Float->Int" floor = floorFloatInt
+"ceiling/Float->Int" ceiling = ceilingFloatInt
+"round/Float->Int" round = roundFloatInt
+ #-}
instance RealFrac Float where
- {-# SPECIALIZE properFraction :: Float -> (Int, Float) #-}
- {-# SPECIALIZE round :: Float -> Int #-}
-
- {-# SPECIALIZE properFraction :: Float -> (Integer, Float) #-}
- {-# SPECIALIZE round :: Float -> Integer #-}
-
-- ceiling, floor, and truncate are all small
- {-# INLINE ceiling #-}
- {-# INLINE floor #-}
- {-# INLINE truncate #-}
+ {-# INLINE [1] ceiling #-}
+ {-# INLINE [1] floor #-}
+ {-# INLINE [1] truncate #-}
-- We assume that FLT_RADIX is 2 so that we can use more efficient code
#if FLT_RADIX != 2
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
(m,_) -> encodeFloat m (negate (floatDigits x))
scaleFloat k x = case decodeFloat x of
- (m,n) -> encodeFloat m (n+k)
+ (m,n) -> encodeFloat m (n + clamp bf k)
+ where bf = FLT_MAX_EXP - (FLT_MIN_EXP) + 4*FLT_MANT_DIG
+
isNaN x = 0 /= isFloatNaN x
isInfinite x = 0 /= isFloatInfinite x
isDenormalized x = 0 /= isFloatDenormalized x
instance Show Float where
showsPrec x = showSignedFloat showFloat x
- showList = showList__ (showsPrec 0)
+ showList = showList__ (showsPrec 0)
\end{code}
%*********************************************************
instance Real Double where
- toRational x = (m%1)*(b%1)^^n
- where (m,n) = decodeFloat x
- b = floatRadix x
+ toRational (D# x#) =
+ case decodeDoubleInteger x# of
+ (# m, e# #)
+ | e# >=# 0# ->
+ shiftLInteger m e# :% 1
+ | (int2Word# (toInt# 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 x = fromRat x
+ 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
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 #-}
+-- RULES for Integer and Int
+{-# RULES
+"properFraction/Double->Integer" properFraction = properFractionDoubleInteger
+"truncate/Double->Integer" truncate = truncateDoubleInteger
+"floor/Double->Integer" floor = floorDoubleInteger
+"ceiling/Double->Integer" ceiling = ceilingDoubleInteger
+"round/Double->Integer" round = roundDoubleInteger
+"properFraction/Double->Int" properFraction = properFractionDoubleInt
+"truncate/Double->Int" truncate = double2Int
+"floor/Double->Int" floor = floorDoubleInt
+"ceiling/Double->Int" ceiling = ceilingDoubleInt
+"round/Double->Int" round = roundDoubleInt
+ #-}
instance RealFrac Double where
- {-# SPECIALIZE properFraction :: Double -> (Int, Double) #-}
- {-# SPECIALIZE round :: Double -> Int #-}
-
- {-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
- {-# SPECIALIZE round :: Double -> Integer #-}
-
-- ceiling, floor, and truncate are all small
- {-# INLINE ceiling #-}
- {-# INLINE floor #-}
- {-# INLINE truncate #-}
+ {-# INLINE [1] ceiling #-}
+ {-# INLINE [1] floor #-}
+ {-# INLINE [1] 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)
+ (fromInteger m * 2 ^ n, 0.0)
else
- case (quotRem m (b^(negate n))) of { (w,r) ->
+ case (quotRem m (2^(negate n))) of { (w,r) ->
(fromInteger w, encodeFloat r n)
}
}
(m,_) -> encodeFloat m (negate (floatDigits x))
scaleFloat k x = case decodeFloat x of
- (m,n) -> encodeFloat m (n+k)
+ (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
instance Show Double where
showsPrec x = showSignedFloat showFloat x
- showList = showList__ (showsPrec 0)
+ showList = showList__ (showsPrec 0)
\end{code}
%*********************************************************
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}
-- | Show a signed 'RealFloat' value to full precision
--- using standard decimal notation for arguments whose absolute value lies
+-- 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)
| 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
+ where
base = 10
doFmt format (is, e) =
-- 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.
+-- and returns a list of digits and an exponent.
-- In particular, if @x>=0@, and
--
-- > floatToDigits base x = ([d1,d2,...,dn], e)
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
- -- f :: Integer, log :: Float -> Float,
+ -- f :: Integer, log :: Float -> Float,
-- ceiling :: Float -> Int
ceiling ((log (fromInteger (f+1) :: Float) +
fromIntegral e * log (fromInteger b)) /
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
fromRat x = x'
where x' = f e
--- If the exponent of the nearest floating-point number to x
+-- 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
\begin{code}
-- | Converts a 'Rational' value into any type in class 'RealFloat'.
-{-# SPECIALISE fromRat :: Rational -> Double,
- Rational -> Float #-}
+{-# RULES
+"fromRat/Float" fromRat = (fromRational :: Rational -> Float)
+"fromRat/Double" fromRat = (fromRational :: Rational -> Double)
+ #-}
fromRat :: (RealFloat a) => Rational -> a
-- Deal with special cases first, delegating the real work to fromRat'
-- 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.
+-- be very slow! We are just slightly more clever, except for base 2, where
+-- we take advantage of the representation of Integers.
+-- The general case could be improved by a lookup table for
+-- approximating the result by integerLog2 i / integerLog2 b.
integerLogBase :: Integer -> Integer -> Int
integerLogBase b i
| i < b = 0
- | otherwise = doDiv (i `div` (b^l)) l
- where
- -- Try squaring the base first to cut down the number of divisions.
- l = 2 * integerLogBase (b*b) i
+ | 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.
- doDiv :: Integer -> Int -> Int
- doDiv x y
- | x < b = y
- | otherwise = doDiv (x `div` b) (y+1)
+The below is an adaption of fromRat' for the conversion to
+Float or Double exploiting the know floatRadix and avoiding
+divisions as much as possible.
+\begin{code}
+{-# SPECIALISE fromRat'' :: Int -> Int -> Integer -> Integer -> Float,
+ Int -> Int -> Integer -> Integer -> Double #-}
+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#) ->
+ if ln# <# md#
+ then encodeFloat (n `shiftL` (I# (md# -# 1# -# ln#)))
+ (I# (ln# +# 1# -# ld# -# md#))
+ 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'# ># (ln# +# 1#) -> encodeFloat 0 0
+ | ld'# ==# (ln# +# 1#) ->
+ case integerLog2IsPowerOf2# n of
+ (# _, 0# #) -> encodeFloat 0 0
+ (# _, _ #) -> encodeFloat 1 (minEx - mantDigs)
+ | ld'# <=# 0# ->
+ encodeFloat n (I# ((me# -# md#) -# ld'#))
+ | otherwise ->
+ 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' (minEx-mantDigs)
+ else encodeFloat (n' + 1) (minEx-mantDigs)
+ _ -> encodeFloat (n' + 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
+ | p <= minEx-mantDigs = (p,a,b)
+ | a < (b `shiftL` (mantDigs-1)) = (p-1, a `shiftL` 1, b)
+ | (b `shiftL` mantDigs) <= a = (p+1, a, b `shiftL` 1)
+ | otherwise = (p, a, b)
+ (p', n'', d'') = scale (p0-mantDigs) 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}
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)
-
-int2Float :: Int -> Float
-int2Float (I# x) = F# (int2Float# x)
-
expFloat, logFloat, sqrtFloat :: Float -> Float
sinFloat, cosFloat, tanFloat :: Float -> Float
asinFloat, acosFloat, atanFloat :: Float -> Float
ltDouble (D# x) (D# y) = x <## y
leDouble (D# x) (D# y) = x <=## y
-double2Int :: Double -> Int
-double2Int (D# x) = I# (double2Int# x)
-
-int2Double :: Int -> Double
-int2Double (I# x) = D# (int2Double# x)
-
double2Float :: Double -> Float
double2Float (D# x) = F# (double2Float# x)
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
+
+A few more details in Trac #2251, and the patch message
"Add RULES for realToFrac from Int".
%*********************************************************
= 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