Haskell FFI / C MPFR library wrapper woes

Question

In order to create an arbitrary precision floating point / drop in replacement for Double, I'm trying to wrap MPFR 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 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.

Questions

• 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 code "works" work? Does the Haskell code "noworks" work?) Can anyone on Linux and Windows attempt to compile/run and see if you get the same results?

C code: (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);
    }  
}

Haskell code: (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]

Answer 1

Judah Jacobsen answered this on the Haskell-cafe mailing list:

This is a known issue with GHC because of the way GHC uses GMP internally (to maintain Integers).

Apparently C data in the heap is left alone by GHC in basically all cases except code that uses the FFI to access GMP or any C library that relies on GMP (like MPFR that I wanted to use). There are some workarounds (painful) but the "right" way would be to either hack GHC (hard) or get the Simons to remove GHC's dependence on GMP (harder).

Answer 2

I see the problem too, on a

$ uname -a
Linux burnup 2.6.26-gentoo-r1 #1 SMP PREEMPT Tue Sep 9 00:05:54 EDT 2008 i686 Intel(R) Pentium(R) 4 CPU 2.80GHz GenuineIntel GNU/Linux
$ gcc --version
gcc (GCC) 4.2.4 (Gentoo 4.2.4 p1.0)
$ ghc --version
The Glorious Glasgow Haskell Compilation System, version 6.8.3

I also see the output changing from 1.0000...000 to 1.0000...[garbage].

Let's see, the following does work:

main = do
    setPrecisionDecimal 640
    mapM_ (const $ c_set_signed_int (fromInteger 1) >>= c_show) [1..10000]

which narrows down the problem to parts of one being somehow clobbered during runtime. Looking at the output of ghc -C and ghc -S, though, isn't giving me any hints.

Hmm, ./noworks +RTS -H1G also works, and ./noworks +RTS -k[n]k, for varying values of [n], demonstrate failures in different ways.

I've got no solid leads, but there are two possibilities that jump to my mind:

• GMP, which the GHC runtime uses, and MPFR having some weird interaction
• stack space for C functions called within the GHC runtime is limited, and MPFR not dealing well

That being said... is there a reason you're rolling your own bindings rather than use HMPFR?

Answer 3

Aleš Bizjak, maintainer of HMPFR posted to haskell-cafe and showed how to keep GHC from controlling allocation of the limbs (and hence leaving them alone, instead of GCing them and clobbering them):

mpfr_ptr mpf_new_mpfr()  
{  
  mpfr_ptr result = malloc(sizeof(__mpfr_struct));  
  if (result == NULL) return NULL;  
  /// these three lines:  
  mp_limb_t * limb = malloc(mpfr_custom_get_size(mpfr_get_default_prec()));  
  mpfr_custom_init(limb, mpfr_get_default_prec());  
  mpfr_custom_init_set(result, MPFR_NAN_KIND, 0, mpfr_get_default_prec(), limb);  
  return result;  
}

To me, this is much easier than joining the effort to write a replacement for GMP in GHC, which would be the only alternative if I really wanted to use any library that depends on GMP.