GMP and conversion of rational to doubles

Kish Shen
2009-03-09
2013-06-06
  • Kish Shen
    Kish Shen
    2009-03-09

    Hi,

    The code I am porting to Win64 using MinGW-W64 uses GMP, so I built GMP 4.2.2 with MinGW-W64, after applying the patch I found in the toolchain patches. I crossed compiled this code on a 32 bit Linux system, so I was unable to run the check (normally done with 'make check') on the build.

    The GMP seems to work with the simple cases I tested it with, but when I ran the test suite for our code on it, I found some problems related to the conversion of GMP rationals to double. For some rationals that cannot be represented precisely with base 2, the double that is returned is not the cloest double to the rational, for example, for the
    rational 4/5, the GMP for W64 returns 0.79999999999999993 rather than 0.80000000000000004, which is a closer double to 0.8 (and is returned by the same GMP on other platforms, with the same double representation).

    Is this simply due to some default rounding mode for Win64? (Our Win32 code, compiled with the 32 bit MinGW, does return 0.80000000000000004, although we are currently using GMP 3.3 there) Is it a known problem and is there a work-around? Alternatively, is there some pre-compiled Win64 GMP that could be used with MinGW-W64 compiled code?

    Thanks in advance for any information/help!

    Kish

     
    • Kai Tietz
      Kai Tietz
      2009-03-09

      Yes, at the moment our code for setting up custom rounding modes via fesetround isn't completely present for amd64. There is still some work necessary. I assume I could fix that this week. Or if you want you can sent a patch for it (see in our crt misc/fesetround.c - we don't set the SSE2 control words for rounding).

      Cheers,
      Kai

       
    • Kish Shen
      Kish Shen
      2009-03-09

      Thanks for your quick reply.

      As I understand it, fesetround is a Unix procedure -- we use it for rounding control on Mac OS X (we do our own rounding control for our double floating point operations), but we use _controlfp for rounding control on Windows. I don't know the GMP source, but I would have expected it to use platform-specific rounding control procedures like us, but am I correct in thinking that GMP need to use fesetround() even on Windows, at least when compiled by MinGW?

      Cheers,

      Kish

       
      • Kai Tietz
        Kai Tietz
        2009-03-10

        Yes, fe... family is ported to mingw from the *nix world.
        I fixed the fesetround method to support SSE, too (see revision 667).

        Cheers,
        Kai

         
    • Kish Shen
      Kish Shen
      2009-03-12

      Hi Kai,

      I downloaded  and rebuilt MinGW-W64 today, and recompiled GMP using it, and it still seem to produce the same rounding issues with rationals.

      Do you expect the change to fesetround you checked in to fix the rounding issue? Did I download the correct version -- the crt and headers were downloaded using svn from the trunk branch.

      Thanks and cheers,

      Kish

       
      • Kai Tietz
        Kai Tietz
        2009-03-14

        Hmm, well. I am not sure if GMP uses fe functions, but the patch fix the rounding via fe functions.

        I'll take a closer look. I am not sure if in all cases fninit is enough. Possibly we have here the troubles, that by default fpu and sse have different settings. (see CRT_fp10.c) and fesetenv.)

        Cheers,
        Kai

         
    • Kish Shen
      Kish Shen
      2009-03-16

      Hi,

      I have done some searching for the rounding issue for rational to double in GMP (mpq_get_d()), and it seems that from GMP 4.2, it is defined as rounding towards zero, so rounding 4/5 to 0.79999999999999993 may be the correct specified behaviour. We have been using GMP 4.1 or earlier in our code, where the rounding behaviour of mpq_get_d() was not defined, and the behaviour was apparently rounding to the nearest double, at least on the platforms we are using.

      So the real question may be why GMP 4.2.2 is behaving differently with 32 bit Linux, where 4/5 gives 0.80000000000000004 -- I have just checked that this is indeed the case. I am planning to check with the GMP mailing list to see what they say.

      Cheers,

      Kish