Floating point manipulation : ulp and coerce IEEE-754 Double# into Word64#

[malosh]

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

[simonmar]

What I had in mind was something like:

coerceDoubleToWord64# :: Double# -> Word64#
coerceFloatToWord32#  :: Float#  -> Word#

(note we don't have a Word32# type)

By the way, you can do this without the FFI right now, by poking the Double into memory, and then peeking it out as a Word64. Not terribly efficient, but better than using the FFI.

[daniel.is.fischer]

I would like to see those coercions and their inverses, too.

[SimonMeier]

I would like to see those coercions and their inverses, too.

[simonmar]

partial patch for this attached. I ran into difficulty because there's no way to allocate temporary memory for the conversion in the native code generator. There are various ways to hack around this (use a spare slot in StgReg?, allocate a new word in the data segment) but none are very nice. We really just want to allocate a slot from the spill area, but that requires cooperation from the register allocator.