Generic get_d_2exp failures

[Niels Möller]

Marc Glisse <marc.glisse at inria.fr> writes:

> 2) The rounding occurs in the addition in
>
>       weight = 1/MP_BASE_AS_DOUBLE;
>       d = up[size - 1];
>       for (i = size - 2; i >= 0; i--)
>         {
>           d += up[i] * weight;
>           weight /= MP_BASE_AS_DOUBLE;
>           if (weight == 0)
>             break;
>         }
>
> we could make that code uglier to make sure there is no rounding up,
> maybe comparing d+up[i]*weight to d+oldweight, and trying again by
> zeroing out an increasing number of low bits of up[i].

It would be nice if we could find a portable way to add two floating
point values without rounding up. Would something like this work?

  s = a + b; /* Assume a > b */
  r = (s - a) - b;  /* No rounding expected here. */
  if (r > 0)
    s -= 2*r;

But doesn't make the tests pass  in my initial testing. I guess it might
be necessary to isolate the high bit of r (possibly be converting to
mp_limb_t ans usign count_leading_zeros), which should correspond to
the lsb of the mantissa (or if it's half the lsb).

And might also need special handling of the case that s, after rounding,
is a power of two.

The variant I tested:

@@ -368,7 +369,13 @@ mpn_get_d (mp_srcptr up, mp_size_t size,
       d = up[size - 1];
       for (i = size - 2; i >= 0; i--)
 	{
-	  d += up[i] * weight;
+	  double l, s, r;
+	  l = up[i] * weight;
+	  s = d + l;
+	  r = (s - d) - l;
+	  if (r > 0)
+	    s -= 2*r;
+	  d = s;
 	  weight /= MP_BASE_AS_DOUBLE;
 	  if (weight == 0)
 	    break;

> 3) Drop the generic path. It hasn't passed the testsuite for a long
> time, can't be that important.

If we keep it, we'd need to (i) fix it, and (ii) test it regularly.

Regards,
/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.