lhc-0.6.20090126: lib/base/src/Prelude/Float.hs
{-# OPTIONS_LHC -N -fffi -fm4 #-}
module Prelude.Float(readDouble,doubleToDigits,doubleToRational) where
import Lhc.Order
import Lhc.Basics
import Lhc.Monad
import Lhc.IO
import Lhc.Float
import Lhc.Num
import Lhc.Types
import Data.Word
import Foreign.Storable
import Foreign.C.Types
import Foreign.Ptr
import Foreign.Marshal.Alloc
import Numeric
import Prelude.Text
import Lhc.List
import Prelude((^),(^^),elem,take)
m4_define(INST,{{
foreign import primitive "FDiv" divide$2 :: $2 -> $2 -> $2
foreign import primitive "FPwr" exponent$2 :: $2 -> $2 -> $2
foreign import primitive "FAtan2" atan2$1 :: $1 -> $1 -> $1
foreign import primitive "F2I" toInteger$1 :: $1 -> Integer
foreign import primitive "const.M_PI" c_pi$1 :: $1
instance Fractional $1 where
$1 x / $1 y = $1 (divide$2 x y)
fromRational x = fromInteger (numerator x) / fromInteger (denominator x)
m4_define(FI,foreign import primitive "{{$}}1" {{$}}2$2 :: $2 -> $2)
FI(Sqrt,sqrt)
FI(Exp,exp)
FI(Log,log)
FI(Sin,sin)
FI(Cos,cos)
FI(Tan,tan)
FI(Sinh,sinh)
FI(Cosh,cosh)
FI(Tanh,tanh)
FI(Asin,asin)
FI(Acos,acos)
FI(Atan,atan)
m4_undefine({{FI}})
instance Floating $1 where
pi = c_pi$1
sqrt ($1 x) = $1 (sqrt$2 x)
exp ($1 x) = $1 (exp$2 x)
log ($1 x) = $1 (log$2 x)
sin ($1 x) = $1 (sin$2 x)
cos ($1 x) = $1 (cos$2 x)
tan ($1 x) = $1 (tan$2 x)
asin ($1 x) = $1 (asin$2 x)
acos ($1 x) = $1 (acos$2 x)
atan ($1 x) = $1 (atan$2 x)
sinh ($1 x) = $1 (sinh$2 x)
cosh ($1 x) = $1 (cosh$2 x)
tanh ($1 x) = $1 (tanh$2 x)
$1 x ** $1 y = $1 (exponent$2 x y)
asinh = c_asinh$1
acosh = c_acosh$1
atanh = c_atanh$1
instance RealFrac $1 where
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 = fromInteger (toInteger$1 x)
round x = fromInteger (toInteger$1 (roundf x))
ceiling x = fromInteger (toInteger$1 (ceilingf x))
floor x = fromInteger (toInteger$1 (floorf x))
properFractionf x = (c_trunc$1 x,x - c_trunc$1 x)
truncatef x = c_trunc$1 x
roundf x = c_nearbyint$1 x
ceilingf x = c_ceil$1 x
floorf x = c_floor$1 x
foreign import ccall "-lm math.h asinh$3" c_asinh$1 :: $1 -> $1
foreign import ccall "-lm math.h acosh$3" c_acosh$1 :: $1 -> $1
foreign import ccall "-lm math.h atanh$3" c_atanh$1 :: $1 -> $1
foreign import ccall "-lm math.h trunc$3" c_trunc$1 :: $1 -> $1
foreign import ccall "-lm math.h ceil$3" c_ceil$1 :: $1 -> $1
foreign import ccall "-lm math.h floor$3" c_floor$1 :: $1 -> $1
foreign import ccall "-lm math.h nearbyint$3" c_nearbyint$1 :: $1 -> $1
foreign import ccall "math.h isnan" c_isnan$3 :: $1 -> CInt
foreign import ccall "math.h isinf" c_isinfinite$3 :: $1 -> CInt
foreign import ccall "math.h signbit" c_signbit$3 :: $1 -> CInt
foreign import ccall "math.h ldexp$3" c_ldexp$3 :: $1 -> CInt -> $1
foreign import ccall "math.h frexp$3" c_frexp$3 :: $1 -> Ptr CInt -> IO $1
}})
INST(Float,Float32_,f)
INST(Double,Float64_)
instance Real Float where
toRational x = (m:%1)*(b:%1)^^n
where (m,n) = decodeFloat x
b = floatRadix x
toDouble x = floatToDouble x
instance Real Double where
toRational x = doubleToRational x
toDouble x = x
instance RealFloat Float where
floatRadix _ = 2
floatDigits _ = 24
floatRange _ = (-125,128)
exponent x = case decodeFloatf x of (_,n) -> n
significand x = case decodeFloatf x of (m,_) -> m
isNaN x = c_isnanf x /= 0
isInfinite x = c_isinfinitef x /= 0
isDenormalized _ = False
isNegativeZero x = x == 0 && c_signbitf x /= 0
isIEEE _ = True
scaleFloat k x = c_ldexpf x (fromInt k)
decodeFloatf x = unsafePerformIO $ alloca $ \ptr -> do
x' <- c_frexpf x ptr
exp <- peek ptr
return (x', fromIntegral exp)
encodeFloat i e = c_ldexpf (fromInteger i) (fromInt e)
decodeFloat x = unsafePerformIO $ alloca $ \ptr -> do
x' <- c_frexp (floatToDouble x) ptr
exp <- peek ptr
let x'' = c_ldexp x' (fromInt $ floatDigits x)
return (double2integer x'', fromIntegral exp - floatDigits x)
atan2 = atan2Float
instance RealFloat Double where
floatRadix _ = 2
floatDigits _ = 53
floatRange _ = (-1021,1024)
exponent x = case decodeFloatf x of (_,n) -> n
significand x = case decodeFloatf x of (m,_) -> m
isNaN x = c_isnan x /= 0
isInfinite x = c_isinfinite x /= 0
isDenormalized _ = False
isNegativeZero x = x == 0 && c_signbit x /= 0
isIEEE _ = True
scaleFloat k x = c_ldexp x (fromInt k)
decodeFloatf x = unsafePerformIO $ alloca $ \ptr -> do
x' <- c_frexp x ptr
exp <- peek ptr
return (x', fromIntegral exp)
encodeFloat i e = c_ldexp (integer2double i) (fromInt e)
decodeFloat x = unsafePerformIO $ alloca $ \ptr -> do
x' <- c_frexp x ptr
exp <- peek ptr
let x'' = c_ldexp x' (fromInt $ floatDigits x)
return (double2integer x'', fromIntegral exp - floatDigits x)
atan2 = atan2Double
foreign import primitive "I2F" integer2float :: Integer -> Float
foreign import primitive "I2F" integer2double :: Integer -> Double
foreign import primitive "F2I" double2integer :: Double -> Integer
readDouble :: ReadS Double
readDouble r = [((fromInteger n * (10^^(k-d))),t) | (n,d,s) <- readFix r,(k,t) <- readExp s] ++
[ (0/0, t) | ("NaN",t) <- lex r] ++
[ (1/0, t) | ("Infinity",t) <- lex r]
where
readFix r = [(read (take 15 $ ds++ds'), length ds', t)
| (ds,d) <- lexDigits r,
(ds',t) <- lexFrac d ]
lexFrac ('.':ds) = lexDigits ds
lexFrac s = [("",s)]
readExp (e:s) | e `elem` "eE" = readExp' s
readExp s = [(0,s)]
readExp' ('-':s) = [(-k,t) | (k,t) <- readDec s]
readExp' ('+':s) = readDec s
readExp' s = readDec s
doubleToDigits :: Integer -> Double -> ([Int], Int)
doubleToDigits n d | n `seq` d `seq` d == 0 = ([], 0)
doubleToDigits base' x =
let (f0', e0) = decodeFloat x
base, f0 :: WordMax
base = fromInteger base'
f0 = fromIntegral f0'
(minExp0, _) = floatRange x
p = floatDigits x
b :: WordMax
b = fromInteger $ floatRadix 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 :: WordMax
e :: Int
(f, e) = let n = minExp - e0
in if n > 0 then (f0 `div` (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)
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)
else
(f*2, b^(-e)*2, 1, 1)
k =
let 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
else
ceiling ((log ((fromIntegral (f+1))::Double) +
fromIntegral e * log (fromIntegral b)) /
log (fromIntegral 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::Int)
gen ds rn sN mUpN mDnN | rn `seq` sN `seq` mUpN `seq` mDnN `seq` True =
let (dn, rn') = (rn * base) `divMod` sN
mUpN' = mUpN * base
mDnN' = mDnN * base
in case (rn' < mDnN', rn' + mUpN' > sN) of
(True, False) -> toInt dn : ds
(False, True) -> toInt (dn+1) : ds
(True, True) -> if rn' * 2 < sN then toInt dn : ds else toInt (dn+1) : ds
(False, False) -> gen (toInt dn:ds) rn' sN mUpN' mDnN'
rds,rrds :: [Int]
rrds = reverse rds
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)
expt :: WordMax -> Int -> WordMax
expt base n = base^n
in k `seq` f `seq` e `seq` b `seq` rrds `seq` (rrds, k)
doubleToRational :: Double -> Rational
doubleToRational x = (m:%1)*(b:%1)^^n where
(m,n) = decodeFloat x
b = floatRadix x