MPFR / FFI - Nefarious (Simple?) bug

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

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
_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: MPFR / FFI - Nefarious (Simple?) bug

Chris Mears
"Jared Updike" <[hidden email]> writes:

> 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?

I can't help you fix the problem, but can confirm that it behaves as you
describe on my Linux system.  The program "works" produces correct
output, but noworks gives:

==================================================
[...]
290
1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
291
1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000315821571076304370386829039486164636536120569099678825575157477354990710405658927062619547595122867857873765133026050519075772712996650488675128428778696175272153632855303857710094714398932056909218681833163883482512398576534211955878126971368731378321605715030861892325510831739618679352919703855747262870175016933079782350322237289103929997053812694251557912158924211585973481445271933347978517620
292
[...]
384
1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000315821571076304370386829039486164636536120569099678825575157477354990710405658927062619547595122867857873765133026050519075772712996650488675128428778696175272153632855303857710094714398932056909218681833163883482512398576534211955878126971368731378321605715030861892325510831739618679352919703855747262870175016933079782350322237289103929997053812694251557912158924211585973481445271933347978517620
385
1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000315821571076304370386829039486164636536120569099678825575157477354990710405658927062619547595122867857873765133026050519075772712996650488675128428778696175272153632855303857710094714398932056909218681833163883482512398576534211955878126971368731378321605715030861892281684860507747503219323326708248193093246881500522519352192987169395899758472958160241015525770811984166890239749344848735134541218
386
../get_str.c:149:  assertion failed: size_s1 >= m
Aborted
==================================================

All the values from 291-384 are the same, but 385 is slightly different.
_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: MPFR / FFI - Nefarious (Simple?) bug

Alexander Dunlap
On Sun, Oct 5, 2008 at 5:38 PM, Chris Mears <[hidden email]> wrote:

> "Jared Updike" <[hidden email]> writes:
>
>> 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?
>
> I can't help you fix the problem, but can confirm that it behaves as you
> describe on my Linux system.  The program "works" produces correct
> output, but noworks gives:
>
> ==================================================
> [...]
> 290
> 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
> 291
> 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000315821571076304370386829039486164636536120569099678825575157477354990710405658927062619547595122867857873765133026050519075772712996650488675128428778696175272153632855303857710094714398932056909218681833163883482512398576534211955878126971368731378321605715030861892325510831739618679352919703855747262870175016933079782350322237289103929997053812694251557912158924211585973481445271933347978517620
> 292
> [...]
> 384
> 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000315821571076304370386829039486164636536120569099678825575157477354990710405658927062619547595122867857873765133026050519075772712996650488675128428778696175272153632855303857710094714398932056909218681833163883482512398576534211955878126971368731378321605715030861892325510831739618679352919703855747262870175016933079782350322237289103929997053812694251557912158924211585973481445271933347978517620
> 385
> 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000315821571076304370386829039486164636536120569099678825575157477354990710405658927062619547595122867857873765133026050519075772712996650488675128428778696175272153632855303857710094714398932056909218681833163883482512398576534211955878126971368731378321605715030861892281684860507747503219323326708248193093246881500522519352192987169395899758472958160241015525770811984166890239749344848735134541218
> 386
> ../get_str.c:149:  assertion failed: size_s1 >= m
> Aborted
> ==================================================
>
> All the values from 291-384 are the same, but 385 is slightly different.
> _______________________________________________
> Haskell-Cafe mailing list
> [hidden email]
> http://www.haskell.org/mailman/listinfo/haskell-cafe
>

My results (from Linux i686) are slightly different. I am getting

[...]
287
1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
288
get_str.c:149:  assertion failed: size_s1 >= m
Aborted

i.e. I never get the random numbers; they are all zeroes but it aborts
after 287 successful repetitions.

Changing the default heap size (with +RTS -H<size>) changes how far it
runs (running with -H32M let it go past 8000 repetitions), so this
might be a garbage-collection issue.

Hope that is some help to you.

Alex
_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: MPFR / FFI - Nefarious (Simple?) bug

Jared Updike
Thanks for trying out my code. I'm glad it's not particular to my system.

I suspected it had to do with the GHC RTS monkeying around with the
heap (lower precisions print more iterations and maybe aren't moving
through the heap as fast?) but I think when I added a statement to
print out the hex address for each Ptr MPFR that those pointers (to
'one') in particular weren't getting moved. Likely the stuff living in
the heap, pointed at in the MPFR struct, specifically the pointer to
the limbs (the mantissa or actual digits of the number) are getting
kludged.. I don't really understand why this is the case when I don't
think I'm allocating a lot of the heap or anything, just printing the
same number over and over... Is there a space leak somewhere? (I can
also verify that mpf_new_mpfr is only getting called once per program
run by appending to a file each time it gets called). How do I tell
GHC "leave this mess and anything in the C heap completely alone?"

  Jared.

On 10/5/08, Alexander Dunlap <[hidden email]> wrote:

> On Sun, Oct 5, 2008 at 5:38 PM, Chris Mears <[hidden email]> wrote:
>  > "Jared Updike" <[hidden email]> writes:
>  >
>  >> 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?
>  >
>  > I can't help you fix the problem, but can confirm that it behaves as you
>  > describe on my Linux system.  The program "works" produces correct
>  > output, but noworks gives:
>  >
>  > ==================================================
>  > [...]
>  > 290
>  > 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
>  > 291
>  > 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000315821571076304370386829039486164636536120569099678825575157477354990710405658927062619547595122867857873765133026050519075772712996650488675128428778696175272153632855303857710094714398932056909218681833163883482512398576534211955878126971368731378321605715030861892325510831739618679352919703855747262870175016933079782350322237289103929997053812694251557912158924211585973481445271933347978517620
>  > 292
>  > [...]
>  > 384
>  > 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000315821571076304370386829039486164636536120569099678825575157477354990710405658927062619547595122867857873765133026050519075772712996650488675128428778696175272153632855303857710094714398932056909218681833163883482512398576534211955878126971368731378321605715030861892325510831739618679352919703855747262870175016933079782350322237289103929997053812694251557912158924211585973481445271933347978517620
>  > 385
>  > 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000315821571076304370386829039486164636536120569099678825575157477354990710405658927062619547595122867857873765133026050519075772712996650488675128428778696175272153632855303857710094714398932056909218681833163883482512398576534211955878126971368731378321605715030861892281684860507747503219323326708248193093246881500522519352192987169395899758472958160241015525770811984166890239749344848735134541218
>  > 386
>  > ../get_str.c:149:  assertion failed: size_s1 >= m
>  > Aborted
>  > ==================================================
>  >
>  > All the values from 291-384 are the same, but 385 is slightly different.
>
> > _______________________________________________
>  > Haskell-Cafe mailing list
>  > [hidden email]
>  > http://www.haskell.org/mailman/listinfo/haskell-cafe
>  >
>
>  My results (from Linux i686) are slightly different. I am getting
>
>  [...]
>  287
>  1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
>  288
>
> get_str.c:149:  assertion failed: size_s1 >= m
>  Aborted
>
>
> i.e. I never get the random numbers; they are all zeroes but it aborts
>  after 287 successful repetitions.
>
>  Changing the default heap size (with +RTS -H<size>) changes how far it
>  runs (running with -H32M let it go past 8000 repetitions), so this
>  might be a garbage-collection issue.
>
>  Hope that is some help to you.
>
>  Alex
>
_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: MPFR / FFI - Nefarious (Simple?) bug

Judah Jacobson
On Sun, Oct 5, 2008 at 6:01 PM, Jared Updike <[hidden email]> wrote:

> Thanks for trying out my code. I'm glad it's not particular to my system.
>
> I suspected it had to do with the GHC RTS monkeying around with the
> heap (lower precisions print more iterations and maybe aren't moving
> through the heap as fast?) but I think when I added a statement to
> print out the hex address for each Ptr MPFR that those pointers (to
> 'one') in particular weren't getting moved. Likely the stuff living in
> the heap, pointed at in the MPFR struct, specifically the pointer to
> the limbs (the mantissa or actual digits of the number) are getting
> kludged.. I don't really understand why this is the case when I don't
> think I'm allocating a lot of the heap or anything, just printing the
> same number over and over... Is there a space leak somewhere? (I can
> also verify that mpf_new_mpfr is only getting called once per program
> run by appending to a file each time it gets called). How do I tell
> GHC "leave this mess and anything in the C heap completely alone?"
>

Usually, you can expect GHC to leave the C heap alone.  The exception,
unfortunately, is GMP (which is used by MPFR).  See the following
ticket:

http://hackage.haskell.org/trac/ghc/ticket/311

I'm guessing that's the cause of the corruption.  The relevant note
from that ticket:

> It's a known problem that you can't use GMP via the FFI from
> your Haskell code in GHC, or use a C library that internally
> uses GMP.  We should really use a private version of GMP
> with function names that don't overlap, or better still
> replace GMP altogether.

The comment in that ticket from Benedict Huber details some possible
workarounds.

Best,
-Judah
_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: MPFR / FFI - Nefarious (Simple?) bug

Michal Konečný
In reply to this post by Jared Updike
Jared Updike wrote on  6 Oct 02:25:
> In order to create an arbitrary precision floating point / drop in
> replacement for Double, I'm trying to wrap MPFR
> (http://www.mpfr.org/)

Have you looked at:
http://hackage.haskell.org/cgi-bin/hackage-scripts/package/hmpfr

Michal
--
|o| Michal Konecny
|o|    http://www-users.aston.ac.uk/~konecnym/
|o|    office: (+42) (0)121 204 3462
|o| PGP key http://www-users.aston.ac.uk/~konecnym/ki.aow

_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe

signature.asc (196 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: MPFR / FFI - Nefarious (Simple?) bug

Henning Thielemann
In reply to this post by Jared Updike

On Sun, 5 Oct 2008, Jared Updike wrote:

> 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.

Don't forget to add it to
   http://haskell.org/haskellwiki/Libraries_and_tools/Mathematics#Real_and_rational_numbers
when you are ready.

_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: MPFR / FFI - Nefarious (Simple?) bug

Aleš Bizjak
In reply to this post by 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.

I've done an interface to MPFR. See Michal's message. It isn't  
particularly fast
but it works. You could improve/suggest improvements there, if you are  
interested.

> Why does the C version work, but the Haskell version flake out?

I'm more or less guessing here, but it might have something to do with mpfr
calling __gmp_allocate_func in init2.c which is called inside  
mpfr_init_set_si, coupled
with what Judah Jacobson said. This is the only possible problem I can  
think of because
when using custom memory interface and allocating with malloc all the  
problems seem to disappear.

If you substitute your mpfr_set_signed_int for the one below, the noworks  
should become works.

mpfr_ptr mpf_set_signed_int(int x)
{
   mpfr_ptr result = mpf_new_mpfr();
   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);
   if (result == NULL)
     return NULL;
   mpfr_set_si(result, x, GMP_RNDN);
   return result;
}
_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe
Reply | Threaded
Open this post in threaded view
|

Re: MPFR / FFI - Nefarious (Simple?) bug

Jared Updike
>  If you substitute your mpfr_set_signed_int for the one below, the noworks
> should become works.

Yes, this works perfectly. Thank you for saving my project! (Had I
known about HMPFR when I started this project, I would have just used
that instead of rolling my own, but it looks like HMPFR was uploaded
at the end of September though I started rolling my own bindings in
June or July... you know how life goes: work, moving, vacation, wife
starting grad school, etc.)

Last night I was seriously considering dropping work on my current
project and joining the effort to replace GMP in GHC as the only
solution. But now I

Thanks,
  Jared.
_______________________________________________
Haskell-Cafe mailing list
[hidden email]
http://www.haskell.org/mailman/listinfo/haskell-cafe