fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny (v3)

[Paul Zimmermann]

       Dear Adhemerval,

thank you again for your useful comments.

> Does it require any ULP adjustment (at least on the architecture you
> tested it)?

I reset to 0 the j0/float values in sysdeps/x86_64/fpu/libm-test-ulps,
then ran "make regen-ulps", and got:

$ diff ../sysdeps/x86_64/fpu/libm-test-ulps /localdisk/zimmerma/glibc/build/math/NewUlps
1315c1315
< float: 0
---
> float: 8
1321c1321
< float: 0
---
> float: 4
1327c1327
< float: 0
---
> float: 5
1333c1333
< float: 0
---
> float: 5

which were the original values except for j0_towardzero (it was 6).
However the same applies to "master", so this is independent from my change.
Anyway, I decreased the value from 6 to 5.

Below is a v3.

Paul

>From 34d336cf2e8495bba21254157be13e02f52d934a Mon Sep 17 00:00:00 2001
From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Date: Thu, 30 Jul 2020 09:16:36 +0200
Subject: [PATCH] fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is
 tiny (v3)

---
 math/auto-libm-test-in            |  2 ++
 sysdeps/ieee754/flt-32/e_j0f.c    | 17 ++++++++++++++++-
 sysdeps/x86_64/fpu/libm-test-ulps |  2 +-
 3 files changed, 19 insertions(+), 2 deletions(-)

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 4414e54d93..5d488a8711 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -5748,6 +5748,8 @@ j0 0x1p16382
 j0 0x1p16383
 # the next value generates larger error bounds on x86_64 (binary32)
 j0 0x2.602774p+0 xfail-rounding:ibm128-libgcc
+# the next value exercises the flt-32 code path for x >= 2^127
+j0 0x8.2f4ecp+124
 
 j1 -1.0
 j1 0.0
diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index c89b9f2688..91e8de8fe3 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -55,7 +55,22 @@ __ieee754_j0f(float x)
 		    z = -__cosf(x+x);
 		    if ((s*c)<zero) cc = z/ss;
 		    else	    ss = z/cc;
-		}
+		} else {
+                  /* We subtract (exactly) a value x0 such that cos(x0)+sin(x0)
+                     is very near to 0, and use the identity
+                     sin(x-x0) = sin(x)*cos(x0)-cos(x)*sin(x0) to get
+                     sin(x) + cos(x) with extra accuracy. */
+                  float x0 = 0xe.d4108p+124f;
+                  float y = x - x0; /* exact */
+                  /* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
+                  z = __sinf (y);
+                  float eps = 0x1.5f263ep-24f;
+                  /* cos(x0) ~ -sin(x0) + eps */
+                  z += eps * __cosf (x);
+                  /* now z ~ (sin(x)-cos(x))*cos(x0) */
+                  float cosx0 = -0xb.504f3p-4f;
+                  cc = z / cosx0;
+                }
 	/*
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index d71d86d961..2561243fb7 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -1324,7 +1324,7 @@ ldouble: 4
 
 Function: "j0_towardzero":
 double: 5
-float: 6
+float: 5
 float128: 2
 ldouble: 5
 
-- 
2.27.0