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, 50 insertions, 0 deletions

diff --git a/math/auto-libm-test-out-j0 b/math/auto-libm-test-out-j0
index cc1fea074e..85c4148de3 100644
--- a/math/auto-libm-test-out-j0
+++ b/math/auto-libm-test-out-j0
@@ -1359,6 +1359,56 @@ j0 0x2.602774p+0 xfail-rounding:ibm128-libgcc
 = j0 tonearest ibm128 0x2.602774p+0 : 0x3.e83779fe19911fa806cee83ab7p-8 : inexact-ok
 = j0 towardzero ibm128 0x2.602774p+0 : 0x3.e83779fe19911fa806cee83ab7p-8 : xfail:ibm128-libgcc inexact-ok
 = j0 upward ibm128 0x2.602774p+0 : 0x3.e83779fe19911fa806cee83ab8p-8 : xfail:ibm128-libgcc inexact-ok
+j0 0x1.04c39cp+6
+= j0 downward binary32 0x4.130e7p+4 : -0x6.dd65d8p-16 : inexact-ok
+= j0 tonearest binary32 0x4.130e7p+4 : -0x6.dd65dp-16 : inexact-ok
+= j0 towardzero binary32 0x4.130e7p+4 : -0x6.dd65dp-16 : inexact-ok
+= j0 upward binary32 0x4.130e7p+4 : -0x6.dd65dp-16 : inexact-ok
+= j0 downward binary64 0x4.130e7p+4 : -0x6.dd65d0005a19p-16 : inexact-ok
+= j0 tonearest binary64 0x4.130e7p+4 : -0x6.dd65d0005a18cp-16 : inexact-ok
+= j0 towardzero binary64 0x4.130e7p+4 : -0x6.dd65d0005a18cp-16 : inexact-ok
+= j0 upward binary64 0x4.130e7p+4 : -0x6.dd65d0005a18cp-16 : inexact-ok
+= j0 downward intel96 0x4.130e7p+4 : -0x6.dd65d0005a18c04p-16 : inexact-ok
+= j0 tonearest intel96 0x4.130e7p+4 : -0x6.dd65d0005a18c04p-16 : inexact-ok
+= j0 towardzero intel96 0x4.130e7p+4 : -0x6.dd65d0005a18c038p-16 : inexact-ok
+= j0 upward intel96 0x4.130e7p+4 : -0x6.dd65d0005a18c038p-16 : inexact-ok
+= j0 downward m68k96 0x4.130e7p+4 : -0x6.dd65d0005a18c04p-16 : inexact-ok
+= j0 tonearest m68k96 0x4.130e7p+4 : -0x6.dd65d0005a18c04p-16 : inexact-ok
+= j0 towardzero m68k96 0x4.130e7p+4 : -0x6.dd65d0005a18c038p-16 : inexact-ok
+= j0 upward m68k96 0x4.130e7p+4 : -0x6.dd65d0005a18c038p-16 : inexact-ok
+= j0 downward binary128 0x4.130e7p+4 : -0x6.dd65d0005a18c03fc1c61141bb64p-16 : inexact-ok
+= j0 tonearest binary128 0x4.130e7p+4 : -0x6.dd65d0005a18c03fc1c61141bb6p-16 : inexact-ok
+= j0 towardzero binary128 0x4.130e7p+4 : -0x6.dd65d0005a18c03fc1c61141bb6p-16 : inexact-ok
+= j0 upward binary128 0x4.130e7p+4 : -0x6.dd65d0005a18c03fc1c61141bb6p-16 : inexact-ok
+= j0 downward ibm128 0x4.130e7p+4 : -0x6.dd65d0005a18c03fc1c61141bcp-16 : inexact-ok
+= j0 tonearest ibm128 0x4.130e7p+4 : -0x6.dd65d0005a18c03fc1c61141bcp-16 : inexact-ok
+= j0 towardzero ibm128 0x4.130e7p+4 : -0x6.dd65d0005a18c03fc1c61141bap-16 : inexact-ok
+= j0 upward ibm128 0x4.130e7p+4 : -0x6.dd65d0005a18c03fc1c61141bap-16 : inexact-ok
+j0 0x1.4b7066p+7
+= j0 downward binary32 0xa.5b833p+4 : 0xf.80edp-20 : inexact-ok
+= j0 tonearest binary32 0xa.5b833p+4 : 0xf.80edp-20 : inexact-ok
+= j0 towardzero binary32 0xa.5b833p+4 : 0xf.80edp-20 : inexact-ok
+= j0 upward binary32 0xa.5b833p+4 : 0xf.80ed1p-20 : inexact-ok
+= j0 downward binary64 0xa.5b833p+4 : 0xf.80ed0055a8e58p-20 : inexact-ok
+= j0 tonearest binary64 0xa.5b833p+4 : 0xf.80ed0055a8e6p-20 : inexact-ok
+= j0 towardzero binary64 0xa.5b833p+4 : 0xf.80ed0055a8e58p-20 : inexact-ok
+= j0 upward binary64 0xa.5b833p+4 : 0xf.80ed0055a8e6p-20 : inexact-ok
+= j0 downward intel96 0xa.5b833p+4 : 0xf.80ed0055a8e5c88p-20 : inexact-ok
+= j0 tonearest intel96 0xa.5b833p+4 : 0xf.80ed0055a8e5c89p-20 : inexact-ok
+= j0 towardzero intel96 0xa.5b833p+4 : 0xf.80ed0055a8e5c88p-20 : inexact-ok
+= j0 upward intel96 0xa.5b833p+4 : 0xf.80ed0055a8e5c89p-20 : inexact-ok
+= j0 downward m68k96 0xa.5b833p+4 : 0xf.80ed0055a8e5c88p-20 : inexact-ok
+= j0 tonearest m68k96 0xa.5b833p+4 : 0xf.80ed0055a8e5c89p-20 : inexact-ok
+= j0 towardzero m68k96 0xa.5b833p+4 : 0xf.80ed0055a8e5c88p-20 : inexact-ok
+= j0 upward m68k96 0xa.5b833p+4 : 0xf.80ed0055a8e5c89p-20 : inexact-ok
+= j0 downward binary128 0xa.5b833p+4 : 0xf.80ed0055a8e5c88af9c0edfe12ap-20 : inexact-ok
+= j0 tonearest binary128 0xa.5b833p+4 : 0xf.80ed0055a8e5c88af9c0edfe12a8p-20 : inexact-ok
+= j0 towardzero binary128 0xa.5b833p+4 : 0xf.80ed0055a8e5c88af9c0edfe12ap-20 : inexact-ok
+= j0 upward binary128 0xa.5b833p+4 : 0xf.80ed0055a8e5c88af9c0edfe12a8p-20 : inexact-ok
+= j0 downward ibm128 0xa.5b833p+4 : 0xf.80ed0055a8e5c88af9c0edfe1p-20 : inexact-ok
+= j0 tonearest ibm128 0xa.5b833p+4 : 0xf.80ed0055a8e5c88af9c0edfe14p-20 : inexact-ok
+= j0 towardzero ibm128 0xa.5b833p+4 : 0xf.80ed0055a8e5c88af9c0edfe1p-20 : inexact-ok
+= j0 upward ibm128 0xa.5b833p+4 : 0xf.80ed0055a8e5c88af9c0edfe14p-20 : inexact-ok
 j0 0x8.2f4ecp+124
 = j0 downward binary32 0x8.2f4ecp+124 : 0xd.331efp-68 : inexact-ok
 = j0 tonearest binary32 0x8.2f4ecp+124 : 0xd.331efp-68 : inexact-ok