Generic get_d_2exp failures
Marc Glisse
marc.glisse at inria.fr
Sat Oct 28 19:38:43 UTC 2017
On Sat, 28 Oct 2017, Niels Möller wrote:
> 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;
Not sure where the 2*r is coming from. And I am not very optimistic that
there is such a simple formula that works for any rounding mode, though I
could easily be wrong. There is also the question of how portable exactly
it needs to be. C99 nextafter can be helpful.
> 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;
It may be easier to assert that FLT_RADIX==2 and use DBL_MANT_DIG to avoid
any rounding. There is even code in one of the tests to compute
DBL_MANT_DIG at runtime (preferably once, at startup) if it isn't
available / reliable at compile-time.
>> 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.
--
Marc Glisse
More information about the gmp-devel
mailing list