Generic get_d_2exp failures
nisse at lysator.liu.se
Sat Oct 28 18:49:59 UTC 2017
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)
> 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)
> 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.
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.
More information about the gmp-devel