Floating point manipulation : ulp and coerce IEEE-754 Double# into Word64#
There are currently two ways to compute the ulp of double numbers with GHC :
- Calling a C function, which requires to allocate a pointer. This is way too expensive when using interval arithmetic (for instance), that compute two ULPs at each arithmetic operations.
- Programming it by hand in haskell with GHC primitive operations, which requires using unsafeCoerce# : this does not work in GHC 6.12.2.
unsafeCoerce# should work, and there should be a primitive ulp# function in GHC, operating on Doubles at least. By the way, here is my haskell code using C for computing it :
foreign import ccall unsafe "math.h frexp" c_frexp::CDouble->(Ptr CInt)->IO () foreign import ccall unsafe "math.h ldexp" c_ldexp::CDouble->CInt->IO CDouble ulp::Double->Double ulp x=unsafePerformIO $ do expon<-alloca (\e->do c_frexp (realToFrac x) e peek e) (c_ldexp 0.5 $ expon-52) >>= return.realToFrac