[Haskell-cafe] MPFR / FFI - Nefarious (Simple?) bug

[Jared Updike]

In order to create an arbitrary precision floating point / drop in
replacement for Double, I'm trying to wrap MPFR (http://www.mpfr.org/)
using the FFI but despite all my efforts the simplest bit of code
doesn't work. It compiles, it runs, but it crashes mockingly after
pretending to work for a while. A simple C version of the code happily
prints the number "1" to (640 decimal places) a total of 10,000 times.
The Haskell version, when asked to do the same, silently corrupts (?)
the data after only 289 print outs of "1.0000...0000" and after 385
print outs, it causes an assertion failure and bombs. I'm at a loss
for how to proceed in debugging this since it "should work".

The code can be perused at http://hpaste.org/10923 (and attached) and
downloaded at http://www.updike.org/mpfr-broken.tar.gz

I'm using GHC 6.83 on FreeBSD 6 and GHC 6.8.2 on Mac OS X. Note you
will need MPFR (tested with 2.3.2) installed with the correct paths
(change the Makefile) for libs and header files (along with those from
GMP) to successfully compile this.

Why does the C version work, but the Haskell version flake out? What
else am I missing when approaching the FFI? I tried StablePtrs and had
the exact same results.

Can someone else verify if this is a Mac/BSD only problem by compiling
and running my code? (Does the C executable"works" work? Does the
Haskell executable "noworks" not work?) Can anyone on Linux and
Windows attempt to compile/run and see if you can get the same
results?

  Jared.

P.S. I actually have much more code, having wrapped dozens of
functions from the MPFR library and created instances using
unsafePerformIO and foreignPtrs and all that useful high level
Haskell/pure stuff--- all of which basically works perfectly except
that it doesn't work at all since it crashes mysteriously after
pretending to work at first. I had to narrow it down to something
simple and the "Print "1.00000..." 10,000 times" was the best I could
do. If I can get that working, I will then move on to something like
the factorial of 1000 (printing intermediate values), which also works
in C, but not Haskell, and after that, maybe everything will "just
work" and I can use/publish this library!

----------------- works.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <gmp.h>
#include <mpfr.h>
#include "mpfr_ffi.c"

int main()
{
  int i;
  mpfr_ptr one;

  mpf_set_default_prec_decimal(640);

  one = mpf_set_signed_int(1);
  for (i = 0; i < 10000; i++)
    {
      printf("%d\n", i);
      mpf_show(one);
    }
}



---------------- Main.hs --- Doesn't work!

module Main where

import Foreign.Ptr            ( Ptr, FunPtr )
import Foreign.C.Types        ( CInt, CLong, CULong, CDouble )
import Foreign.StablePtr      ( StablePtr )

data MPFR = MPFR

foreign import ccall "mpf_set_default_prec_decimal"
    c_set_default_prec_decimal          :: CInt -> IO ()
setPrecisionDecimal                     :: Integer -> IO ()
setPrecisionDecimal decimal_digits = do
    c_set_default_prec_decimal (fromInteger decimal_digits)

foreign import ccall "mpf_show"
   c_show                               :: Ptr MPFR -> IO ()

foreign import ccall "mpf_set_signed_int"
   c_set_signed_int                     :: CLong -> IO (Ptr MPFR)

showNums k n = do
   print n
   c_show k

main = do
   setPrecisionDecimal 640
   one <- c_set_signed_int (fromInteger 1)
   mapM_ (showNums one) [1..10000]



-------------- mpfr_ffi.c

// a little bit of wrapper code... this may be where the problem is
// but since the (working) C code uses this and has no problems, I don't
// understand why the Haskell code appears to be corrupting the limbs
(the mantissa)

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <gmp.h>
#include <mpfr.h>

void mpf_show(mpfr_ptr x)
{
  mpfr_out_str(stdout, 10, 0, x, GMP_RNDN);
  printf("\n");
}


#define LOG2_10       3.3219280948873626
#define IEEE_BITS     53
void mpf_set_default_prec_decimal(int dec_digits)
{
  double bits = LOG2_10 * dec_digits;
  int ibits = (int)bits;
  double ddec;
  int dec;
  //printf("DEC_DIGITS = %d\n", dec_digits);
  //printf("IBITS      = %d\n", ibits);
  ddec = ibits / LOG2_10;
  dec = (int)ddec;
  //printf("DEC        = %d\n", dec);
  while (dec < dec_digits)
  {
    ibits++;
    //printf("IBITS      = %d\n", ibits);
    ddec = ibits / LOG2_10;
    dec = (int)ddec;
    //printf("DEC        = %d\n", dec);
  }
  mpfr_set_default_prec(ibits);
}

mpfr_ptr mpf_new_mpfr()
{
  return (mpfr_ptr)malloc(sizeof(__mpfr_struct));
}

mpfr_ptr mpf_set_signed_int(int x)
{
  mpfr_ptr result = mpf_new_mpfr();
  if (result == NULL)
    return NULL;
  mpfr_init_set_si(result, x, GMP_RNDN);
  return result;
}



--------------- Makefile


CC      = gcc
HC      = ghc
INCLUDE = -I/home/private/local/include -I/usr/local/include
-I/opt/local/include
LIBS    = -L/opt/local/lib -L/usr/local/lib

all:	works noworks

works:	mpfr_ffi.o works.o
	$(CC) -o works works.o -lgmp -lmpfr $(INCLUDE) $(LIBS)

works.o:	works.c
	$(CC) -c works.c $(INCLUDE)

mpfr_ffi.o: mpfr_ffi.c
	$(CC) -c mpfr_ffi.c $(INCLUDE)

noworks:	Main.hs mpfr_ffi.o
	$(HC) mpfr_ffi.o -fffi -lgmp -lmpfr $(LIBS) $(INCLUDE) --make Main -o noworks

clean:
	rm *.o *.hi works noworks