[PATCH] Fix the inaccuracy of j0f/j1f/y0f/y1f [BZ #14469, #14470, #14471, #14472]

[Adhemerval Zanella]

On 31/03/2021 03:12, Paul Zimmermann wrote:
> For j0f/j1f/y0f/y1f, the largest error for all binary32
> inputs is reduced to at most 9 ulps for all rounding modes.
> 
> The new code is enabled only when there is a cancellation at the very end of
> the j0f/j1f/y0f/y1f computation, or for very large inputs, thus should not
> give any visible slowdown on average.  Two different algorithms are used:
> 
> * around the first 64 zeros of j0/j1/y0/y1, approximation polynomials of
>   degree 3 are used, computed using the Sollya tool (https://www.sollya.org/)
> 
> * for large inputs, an asymptotic formula from [1] is used
> 
> [1] Fast and Accurate Bessel Function Computation,
>     John Harrison, Proceedings of Arith 19, 2009.
> 
> Inputs yielding the new largest errors are added to auto-libm-test-in,
> and ulps are regenerated for various targets (thanks Adhemerval Zanella).
> 
> Tested on x86_64 with --disable-multi-arch and on powerpc64le-linux-gnu.
> ---
>  math/auto-libm-test-in             |  20 +-
>  math/auto-libm-test-out-j0         |  50 +++
>  math/auto-libm-test-out-j1         |  50 +++
>  math/auto-libm-test-out-y0         |  50 +++
>  math/auto-libm-test-out-y1         |  75 +++++
>  sysdeps/aarch64/libm-test-ulps     |  70 ++--
>  sysdeps/ieee754/flt-32/e_j0f.c     | 515 ++++++++++++++++++++++++++---
>  sysdeps/ieee754/flt-32/e_j1f.c     | 512 ++++++++++++++++++++++++++--
>  sysdeps/powerpc/fpu/libm-test-ulps |  62 ++--
>  sysdeps/s390/fpu/libm-test-ulps    |  68 ++--
>  sysdeps/sparc/fpu/libm-test-ulps   |  68 ++--
>  sysdeps/x86_64/fpu/libm-test-ulps  |  76 ++---
>  12 files changed, 1371 insertions(+), 245 deletions(-)

Paul,

This version still misses the reduce_aux.h file.