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, 31 insertions, 31 deletions

diff --git a/sysdeps/powerpc/fpu/libm-test-ulps b/sysdeps/powerpc/fpu/libm-test-ulps
index 173388b92a..3010e10c75 100644
--- a/sysdeps/powerpc/fpu/libm-test-ulps
+++ b/sysdeps/powerpc/fpu/libm-test-ulps
@@ -1310,50 +1310,50 @@ float128: 1
 ldouble: 3
 
 Function: "j0":
-double: 2
-float: 8
+double: 3
+float: 9
 float128: 2
-ldouble: 2
+ldouble: 5
 
 Function: "j0_downward":
-double: 2
-float: 4
+double: 6
+float: 9
 float128: 4
 ldouble: 12
 
 Function: "j0_towardzero":
-double: 5
-float: 6
+double: 7
+float: 9
 float128: 4
 ldouble: 16
 
 Function: "j0_upward":
-double: 4
-float: 5
+double: 9
+float: 8
 float128: 5
-ldouble: 6
+ldouble: 14
 
 Function: "j1":
-double: 2
-float: 8
+double: 4
+float: 9
 float128: 4
-ldouble: 3
+ldouble: 6
 
 Function: "j1_downward":
 double: 3
-float: 5
+float: 8
 float128: 4
 ldouble: 7
 
 Function: "j1_towardzero":
-double: 3
-float: 2
+double: 4
+float: 8
 float128: 4
 ldouble: 7
 
 Function: "j1_upward":
-double: 3
-float: 4
+double: 9
+float: 9
 float128: 3
 ldouble: 6
 
@@ -1706,49 +1706,49 @@ ldouble: 5
 
 Function: "y0":
 double: 2
-float: 6
+float: 8
 float128: 3
-ldouble: 1
+ldouble: 10
 
 Function: "y0_downward":
 double: 3
-float: 4
+float: 8
 float128: 4
 ldouble: 10
 
 Function: "y0_towardzero":
 double: 3
-float: 3
+float: 8
 float128: 3
-ldouble: 8
+ldouble: 9
 
 Function: "y0_upward":
 double: 2
-float: 5
+float: 8
 float128: 3
 ldouble: 9
 
 Function: "y1":
 double: 3
-float: 2
+float: 9
 float128: 2
 ldouble: 2
 
 Function: "y1_downward":
-double: 3
-float: 2
+double: 6
+float: 8
 float128: 4
-ldouble: 7
+ldouble: 11
 
 Function: "y1_towardzero":
 double: 3
-float: 2
+float: 9
 float128: 2
 ldouble: 9
 
 Function: "y1_upward":
-double: 5
-float: 2
+double: 6
+float: 9
 float128: 5
 ldouble: 9