[Haskell-cafe] floating point operations and representation

[Alfonso Acosta]

On Thu, Mar 13, 2008 at 1:35 AM, Jacob Schwartz <quark at bluespec.com> wrote:
> I have two questions about using the Double data type and the
>  operations in the Floating typeclass on a computer that uses IEEE
>  floating point numbers.
>
>  I notice that the Floating class only provides "log" (presumably log
>  base 'e') and "logBase" (which, in the latest source that I see for
>  GHC is defined as "log y / log x").  However, in C, the "math.h"
>  library provides specific "log2" and "log10" functions, for extra
>  precision.  A test on IEEE computers (x86 and x86-64), shows that for
>  a range of 64-bit "double" values, the answers in C do differ (in the
>  last bit) if you use "log2(x)" and "log10(x)" versus "log (x) /
>  log(2)" and "log(x) / log(10)".
>
>  I am under the restriction that I need to write Haskell programs using
>  Double which mimic existing C/C++ programs or generated data sets, and
>  get the same answers.  (It's silly, but take it as a given
>  requirement.)  If the C programs are using "log2", then I need "log2"
>  in the Haskell, or else I run the risk of not producing the same
>  answers.  My first thought is to import "log2" and "log10" through the
>  FFI.  I was wondering if anyone on Haskell-Cafe has already done this
>  and/or has a better suggestion about how to get specialized "log2" and
>  "log10" (among the many specialized functions that the "math.h"
>  library provides, for better precision -- for now, I'm just concerned
>  with "log2" and "log10").

The  RULES pragma + FFI solution (as suggested by Henning) seems to be
the most sensible one so far.


>  My second question is how to get at the IEEE bit representation for a
>  Double.  I am already checking "isIEEE n" in my source code (and
>  "floatRadix n == 2").  So I know that I am operating on hardware that
>  implements floating point numbers by the IEEE standard.  I would like
>  to get at the 64 bits of a Double.  Again, I can convert to a CDouble
>  and use the FFI to wrap a C function which casts the "double" to a
>  64-bit number and returns it.  But I'm wondering if there's not a
>  better way to do this natively in Haskell/GHC (perhaps some crazy use
>  of the Storable typeclass?).


Posted in a a previous haskell-cafe message [1]

{-# LANGUAGE MagicHash #-}
import Data.Int
import GHC.Exts

foo :: Double -> Int64 --ghc only
foo (D# x) = I64# (unsafeCoerce# x)

[1] http://www.haskell.org/pipermail/haskell-cafe/2008-February/040000.html