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-j1 b/math/auto-libm-test-out-j1
index 6bc3bbe5a4..f808801976 100644
--- a/math/auto-libm-test-out-j1
+++ b/math/auto-libm-test-out-j1
@@ -993,6 +993,56 @@ j1 0x3.ae4b2p+0 xfail-rounding:ibm128-libgcc
 = j1 tonearest ibm128 0x3.ae4b2p+0 : 0xf.d085c66e86f30267f22d6f787cp-8 : inexact-ok
 = j1 towardzero ibm128 0x3.ae4b2p+0 : 0xf.d085c66e86f30267f22d6f787cp-8 : xfail:ibm128-libgcc inexact-ok
 = j1 upward ibm128 0x3.ae4b2p+0 : 0xf.d085c66e86f30267f22d6f788p-8 : xfail:ibm128-libgcc inexact-ok
+j1 0x1.2f28eap+7 xfail-rounding:binary64 xfail-rounding:binary128 xfail-rounding:intel96 xfail-rounding:ibm128-libgcc
+= j1 downward binary32 0x9.79475p+4 : 0x2.49a6fp-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 tonearest binary32 0x9.79475p+4 : 0x2.49a6f4p-16 : inexact-ok
+= j1 towardzero binary32 0x9.79475p+4 : 0x2.49a6fp-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 upward binary32 0x9.79475p+4 : 0x2.49a6f4p-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 downward binary64 0x9.79475p+4 : 0x2.49a6f3fb00346p-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 tonearest binary64 0x9.79475p+4 : 0x2.49a6f3fb00346p-16 : inexact-ok
+= j1 towardzero binary64 0x9.79475p+4 : 0x2.49a6f3fb00346p-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 upward binary64 0x9.79475p+4 : 0x2.49a6f3fb00348p-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 downward intel96 0x9.79475p+4 : 0x2.49a6f3fb003462ap-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 tonearest intel96 0x9.79475p+4 : 0x2.49a6f3fb003462ap-16 : inexact-ok
+= j1 towardzero intel96 0x9.79475p+4 : 0x2.49a6f3fb003462ap-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 upward intel96 0x9.79475p+4 : 0x2.49a6f3fb003462a4p-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 downward m68k96 0x9.79475p+4 : 0x2.49a6f3fb003462ap-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 tonearest m68k96 0x9.79475p+4 : 0x2.49a6f3fb003462ap-16 : inexact-ok
+= j1 towardzero m68k96 0x9.79475p+4 : 0x2.49a6f3fb003462ap-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 upward m68k96 0x9.79475p+4 : 0x2.49a6f3fb003462a4p-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 downward binary128 0x9.79475p+4 : 0x2.49a6f3fb003462a1f15135cec05ap-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 tonearest binary128 0x9.79475p+4 : 0x2.49a6f3fb003462a1f15135cec05cp-16 : inexact-ok
+= j1 towardzero binary128 0x9.79475p+4 : 0x2.49a6f3fb003462a1f15135cec05ap-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 upward binary128 0x9.79475p+4 : 0x2.49a6f3fb003462a1f15135cec05cp-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 downward ibm128 0x9.79475p+4 : 0x2.49a6f3fb003462a1f15135cecp-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 tonearest ibm128 0x9.79475p+4 : 0x2.49a6f3fb003462a1f15135cecp-16 : inexact-ok
+= j1 towardzero ibm128 0x9.79475p+4 : 0x2.49a6f3fb003462a1f15135cecp-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 upward ibm128 0x9.79475p+4 : 0x2.49a6f3fb003462a1f15135cec1p-16 : xfail:binary64 xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+j1 0x1.a1d20ap+6 xfail-rounding:binary128 xfail-rounding:intel96 xfail-rounding:ibm128-libgcc
+= j1 downward binary32 0x6.874828p+4 : -0x3.d6f0b8p-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 tonearest binary32 0x6.874828p+4 : -0x3.d6f0b4p-16 : inexact-ok
+= j1 towardzero binary32 0x6.874828p+4 : -0x3.d6f0b4p-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 upward binary32 0x6.874828p+4 : -0x3.d6f0b4p-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 downward binary64 0x6.874828p+4 : -0x3.d6f0b408112dp-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 tonearest binary64 0x6.874828p+4 : -0x3.d6f0b408112dp-16 : inexact-ok
+= j1 towardzero binary64 0x6.874828p+4 : -0x3.d6f0b408112cep-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 upward binary64 0x6.874828p+4 : -0x3.d6f0b408112cep-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 downward intel96 0x6.874828p+4 : -0x3.d6f0b408112cf3dp-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 tonearest intel96 0x6.874828p+4 : -0x3.d6f0b408112cf3ccp-16 : inexact-ok
+= j1 towardzero intel96 0x6.874828p+4 : -0x3.d6f0b408112cf3ccp-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 upward intel96 0x6.874828p+4 : -0x3.d6f0b408112cf3ccp-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 downward m68k96 0x6.874828p+4 : -0x3.d6f0b408112cf3dp-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 tonearest m68k96 0x6.874828p+4 : -0x3.d6f0b408112cf3ccp-16 : inexact-ok
+= j1 towardzero m68k96 0x6.874828p+4 : -0x3.d6f0b408112cf3ccp-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 upward m68k96 0x6.874828p+4 : -0x3.d6f0b408112cf3ccp-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 downward binary128 0x6.874828p+4 : -0x3.d6f0b408112cf3cdb53ca3ca7088p-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 tonearest binary128 0x6.874828p+4 : -0x3.d6f0b408112cf3cdb53ca3ca7086p-16 : inexact-ok
+= j1 towardzero binary128 0x6.874828p+4 : -0x3.d6f0b408112cf3cdb53ca3ca7086p-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 upward binary128 0x6.874828p+4 : -0x3.d6f0b408112cf3cdb53ca3ca7086p-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 downward ibm128 0x6.874828p+4 : -0x3.d6f0b408112cf3cdb53ca3ca71p-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 tonearest ibm128 0x6.874828p+4 : -0x3.d6f0b408112cf3cdb53ca3ca71p-16 : inexact-ok
+= j1 towardzero ibm128 0x6.874828p+4 : -0x3.d6f0b408112cf3cdb53ca3ca7p-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
+= j1 upward ibm128 0x6.874828p+4 : -0x3.d6f0b408112cf3cdb53ca3ca7p-16 : xfail:binary128 xfail:intel96 xfail:ibm128-libgcc inexact-ok
 j1 min
 = j1 downward binary32 0x4p-128 : 0x1.fffff8p-128 : inexact-ok underflow errno-erange-ok
 = j1 tonearest binary32 0x4p-128 : 0x2p-128 : inexact-ok underflow errno-erange-ok