mini-gmp mpz_{get,set}_d not fully compatible with GMP
Marco Bodrato
bodrato at mail.dm.unipi.it
Tue Mar 6 18:40:11 UTC 2018
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/
More information about the gmp-devel
mailing list