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

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.
Reviewed-by: Adhemerval Zanella  <adhemerval.zanella@linaro.org>

1 files changed, 35 insertions, 35 deletions

diff --git a/sysdeps/aarch64/libm-test-ulps b/sysdeps/aarch64/libm-test-ulps
index 21ff7dc5ce..a0329dea63 100644
--- a/sysdeps/aarch64/libm-test-ulps
+++ b/sysdeps/aarch64/libm-test-ulps
@@ -1066,44 +1066,44 @@ double: 1
 ldouble: 1
 
 Function: "j0":
-double: 2
-float: 8
+double: 3
+float: 9
 ldouble: 2
 
 Function: "j0_downward":
-double: 2
-float: 4
-ldouble: 4
+double: 6
+float: 9
+ldouble: 9
 
 Function: "j0_towardzero":
-double: 5
-float: 6
-ldouble: 4
+double: 7
+float: 9
+ldouble: 9
 
 Function: "j0_upward":
-double: 4
-float: 5
-ldouble: 5
+double: 9
+float: 8
+ldouble: 7
 
 Function: "j1":
-double: 2
-float: 8
+double: 4
+float: 9
 ldouble: 4
 
 Function: "j1_downward":
 double: 3
-float: 5
-ldouble: 4
+float: 8
+ldouble: 6
 
 Function: "j1_towardzero":
-double: 3
-float: 2
-ldouble: 4
+double: 4
+float: 8
+ldouble: 9
 
 Function: "j1_upward":
-double: 3
-float: 4
-ldouble: 3
+double: 9
+float: 9
+ldouble: 9
 
 Function: "jn":
 double: 4
@@ -1364,42 +1364,42 @@ ldouble: 4
 
 Function: "y0":
 double: 2
-float: 6
+float: 8
 ldouble: 3
 
 Function: "y0_downward":
 double: 3
-float: 4
-ldouble: 4
+float: 8
+ldouble: 7
 
 Function: "y0_towardzero":
 double: 3
-float: 3
+float: 8
 ldouble: 3
 
 Function: "y0_upward":
 double: 2
-float: 5
-ldouble: 3
+float: 8
+ldouble: 4
 
 Function: "y1":
 double: 3
-float: 2
-ldouble: 2
+float: 9
+ldouble: 5
 
 Function: "y1_downward":
-double: 3
-float: 2
-ldouble: 4
+double: 6
+float: 8
+ldouble: 5
 
 Function: "y1_towardzero":
 double: 3
-float: 2
+float: 9
 ldouble: 2
 
 Function: "y1_upward":
-double: 5
-float: 2
+double: 6
+float: 9
 ldouble: 5
 
 Function: "yn":