mini-gmp mpz_{get,set}_d not fully compatible with GMP

[Marco Bodrato]

Ciao,

I was trying to write case tests for the (proposed) mini-mpq layer, and I
fell into a possible problem in che mini-gmp mpz layer,when dealing with
mpz<->double conversion.

I tested the following program:

#include <gmp.h>

int
main (int argc, char **argv)
{
  mpz_t  za, zb;
  int i, cp, cm;
  mpz_init (za);
  mpz_init (zb);

  for (i = 0; ++i < 200;)
    {
      volatile double dp, dm;

      mpz_set_ui (za, 0);
      mpz_setbit (za, i);
      mpz_add_ui (za, za, 1); /* za = 2^i + 1 */
      dp = mpz_get_d (za);
      mpz_set_d (zb, dp);
      cp = mpz_cmp (za, zb);

      mpz_sub_ui (za, za, 2); /* za = 2^i - 1 */
      dm = mpz_get_d (za);
      mpz_set_d (zb, dm);
      cm = mpz_cmp (za, zb);
      printf ("i = %i ) cm = %i, cp = %i, dp - dm = %g\n",
               i, cm, cp, dp-dm);
    }
  mpz_clear (za);
  mpz_clear (zb);
  return 0;
}

With GMP I get:

...
i = 51 ) cm = 0, cp = 0, dp - dm = 2
i = 52 ) cm = 0, cp = 0, dp - dm = 2
i = 53 ) cm = 0, cp = 1, dp - dm = 1
i = 54 ) cm = 1, cp = 1, dp - dm = 2
i = 55 ) cm = 1, cp = 1, dp - dm = 4
i = 56 ) cm = 1, cp = 1, dp - dm = 8
i = 57 ) cm = 1, cp = 1, dp - dm = 16
...

I read this results as: (GMP's) mpz_get_d always returns a double that is
rounded towards the same direction (documentation says "rounding towards
zero"); when 2^i gets larger than the precision of a double, 2^i+1 and
2^i-i are traslated to distant values. Distance increases when increasing
the exponent i.

With mini-gmp I get:

...
i = 51 ) cm = 0, cp = 0, dp - dm = 2
i = 52 ) cm = 0, cp = 0, dp - dm = 2
i = 53 ) cm = 0, cp = 1, dp - dm = 1
i = 54 ) cm = -1, cp = 1, dp - dm = 0
i = 55 ) cm = -1, cp = 1, dp - dm = 0
i = 56 ) cm = -1, cp = 1, dp - dm = 0
i = 57 ) cm = -1, cp = 1, dp - dm = 0
...

I read this results as: (mini-gmp's) mpz_get_d returns the nearest double;
when 2^i gets larger than the precision of a double, 2^i+1 and 2^i-i are
traslated the same value (2^i, I suppose).

Maybe this behaviour of mini-gmp changes with the different
implementations of the type double on different architectures? (maybe it's
a compiler bug on my side? But we do not test it in mini-gmp/tests, and it
passes undetected).


I believe this means that the promise "The implemented functions are fully
compatible with the corresponding GMP functions", that we wrote in
mini-gmp/README is not completely satisfied for those functions.

Which solutions do you suggest? I see 3 possible ways:
 - removing the non "fully compatible" functions;
 - adding them to the "with a few exceptions:" section;
 - correct the behaviour of the functions...

Opinions? Other options?

Ĝis,
m

-- 
http://bodrato.it/papers/