Re: [PATCH 1/4] Add a routine for accurate reduction mod (2pi), for j0f/j1f/y0f/y1f.

[Wilco Dijkstra <Wilco.Dijkstra@arm.com>]

Hi Paul,

> reduce_fast is clearly not accurate enough: for x=0x1.96d666p+7, it gives
> y=0x1.947b0da22168cp+8 (and n=-128), and we clearly see that y is not in
> [-PI/4,PI/4] as advertized. What happens is that in the #else part of
> the code that is executed on my machine, r is larger than 2^31, thus the
> cast (int32_t) r overflows, the value of n is wrong, and thus the returned
> value x - n * p->hpi. Maybe a comment should be added about this limitation
> of reduce_fast?

Well the comment says:

/* Fast range reduction using single multiply-subtract.  Return the modulo of
   X as a value between -PI/4 and PI/4 and store the quadrant in NP.
   The values for PI/2 and 2/PI are accessed via P.  Since PI/2 as a double
   is accurate to 55 bits and the worst-case cancellation happens at 6 * PI/4,
   the result is accurate for |X| <= 120.0.  */

Ie. it won't be accurate enough over 120.0 (and yes, it will overflow at around
201 but that's well outside the allowed range). Note also reduce_large only
works for values larger than 2.0. The idea is that they complement each other,
allowing the common case to use the faster reduction - see eg. sinf:

  else if (__glibc_likely (abstop12 (y) < abstop12 (120.0f)))
      x = reduce_fast (x, p, &n);
  else if (abstop12 (y) < abstop12 (INFINITY))
      x = reduce_large (xi, &n);

We could add an inline helper that has this logic and the sign handling so
that it works as a generic range reducer across the full float range.
As Carlos points out, it's not essential for your patch, so it could be done
in another path.

There is a generic range reducer for double precision, however it is very
complex and slow, so it needs to be rewritten to use the reduce_fast algorithm.

> As for reduce_large, I need to run again the exhaustive tests for all four
> functions (j0, j1, y0, y1) and all four rounding modes to see if it is accurate
> enough.

It should be since it gives you x mod PI/4 with absolute error of about 2^-62.

Cheers,
Wilco