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

[Marco Bodrato]

Ciao,

Il Mar, 6 Marzo 2018 11:58 pm, Marc Glisse ha scritto:
> It is also possible to determine it at runtime (comparing 2^n+{0,1,2} for
> instance should work whether we round up, down or to the nearest;
> tests/misc.c compares 2^n+1-2^n to 1, which also looks good), doing it

Well, stealing the tests_dbl_mant_bits code from tests-misc.c, I wrote the
following code for mini-gmp. It seams to work...

Comments?

-----8<----------8<----
static unsigned int
gmp_tests_dbl_mant_bits (void)
{
  unsigned int n;
  volatile double x, y, d;

  n = 1;
  x = d = 1.0;
  do {
    x *= 2.0;
    y = x + d;
    d = y - x;
  } while (d == 1.0 && ++n < sizeof(double) * CHAR_BIT);

  return n;
}

double
mpz_get_d (const mpz_t u)
{
  static unsigned int c = 0;
  int m;
  mp_limb_t l;
  mp_size_t un;
  double x;
  double B = 2.0 * (double) GMP_LIMB_HIGHBIT;

  un = GMP_ABS (u->_mp_size);

  if (un == 0)
    return 0.0;

  if (c == 0)
    c = gmp_tests_dbl_mant_bits ();

  l = u->_mp_d[--un];
  gmp_clz (m, l);
  m = m + c - GMP_LIMB_BITS;
  if (m < 0)
    l &= GMP_LIMB_MAX << -m;

  for (x = l; --un >= 0;)
    {
      x = B*x;
      if (m > 0) {
	l = u->_mp_d[un];
	m -= GMP_LIMB_BITS;
	if (m < 0)
	  l &= GMP_LIMB_MAX << -m;
	x += l;
      }
    }

  if (u->_mp_size < 0)
    x = -x;

  return x;
}
-----8<----------8<----


Ĝis,
m

-- 
http://bodrato.it/