Generic get_d_2exp failures

[Marc Glisse]

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