fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x) is tiny
Paul Zimmermann
Paul.Zimmermann@inria.fr
Tue Jul 28 08:23:05 GMT 2020
> In new code we generally prefer to use hex float constants in such cases
> where a specific floating-point value is wanted.
thank you Joseph. Here is a new version. The maximal error for x >= 2^127
is now 4 ulps (attained for x=1.740713465e+38).
Total: errors=4220511 (0.10%) errors2=393216 maxerr=4 ulp(s)
Paul
>From 6b731f36b1a5badf4704645d0dda40957cedd0db Mon Sep 17 00:00:00 2001
From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Date: Mon, 27 Jul 2020 19:01:18 +0200
Subject: [PATCH 1/2] fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x)
is tiny
---
sysdeps/ieee754/flt-32/e_j0f.c | 16 ++++++++++++++++
1 file changed, 16 insertions(+)
diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index c89b9f2688..f85d8a59e0 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -56,6 +56,22 @@ __ieee754_j0f(float 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 from 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 = 3.153646966e+38f;
+ float y = x - x0; /* exact */
+ /* sin(y) = sin(x)*cos(x0)-cos(x)*sin(x0) */
+ z = __sinf (y);
+ float eps = 8.17583368e-8f;
+ /* cos(x0) ~ -sin(x0) + eps */
+ z += eps * __cosf (x);
+ /* now z ~ (sin(x)-cos(x))*cos(x0) */
+ float cosx0 = -0.707106740f;
+ 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)
--
2.27.0
>From 44124c42fe519c7dcac829160181ba0bb6c8751c Mon Sep 17 00:00:00 2001
From: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Date: Tue, 28 Jul 2020 10:05:38 +0200
Subject: [PATCH 2/2] use hex float constants as advised by Joseph Myers
---
sysdeps/ieee754/flt-32/e_j0f.c | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index f85d8a59e0..883a1a4d13 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -61,15 +61,15 @@ __ieee754_j0f(float x)
is very near from 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 = 3.153646966e+38f;
+ 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 = 8.17583368e-8f;
+ float eps = 0x1.5f263ep-24f;
/* cos(x0) ~ -sin(x0) + eps */
z += eps * __cosf (x);
/* now z ~ (sin(x)-cos(x))*cos(x0) */
- float cosx0 = -0.707106740f;
+ float cosx0 = -0xb.504f3p-4f;
cc = z / cosx0;
}
/*
--
2.27.0
More information about the Libc-alpha
mailing list