about summary refs log tree commit diff
path: root/sysdeps
diff options
context:
space:
mode:
authorPaul Zimmermann <Paul.Zimmermann@inria.fr>2020-08-07 16:14:53 -0300
committerAdhemerval Zanella <adhemerval.zanella@linaro.org>2020-08-07 16:33:13 -0300
commitb7dd366dbeeff3016d0554c08a5f5c8cfe4c31d1 (patch)
treeecba427036d83b41e0b4fcec49f0120d4a424366 /sysdeps
parent1cfb4715288845ebc55ad664421b48b32de9599c (diff)
downloadglibc-b7dd366dbeeff3016d0554c08a5f5c8cfe4c31d1.tar.gz
glibc-b7dd366dbeeff3016d0554c08a5f5c8cfe4c31d1.tar.xz
glibc-b7dd366dbeeff3016d0554c08a5f5c8cfe4c31d1.zip
math: Fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny
Checked on x86_64-linux-gnu and i686-linux-gnu.
Diffstat (limited to 'sysdeps')
-rw-r--r--sysdeps/ieee754/flt-32/e_j0f.c17
1 files changed, 16 insertions, 1 deletions
diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index c89b9f2688..5d29611eb7 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)