From mboxrd@z Thu Jan  1 00:00:00 1970
Return-Path: <azanella@sourceware.org>
Received: by sourceware.org (Postfix, from userid 1791)
 id 26B703857C55; Fri,  7 Aug 2020 19:33:34 +0000 (GMT)
DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org 26B703857C55
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
From: Adhemerval Zanella <azanella@sourceware.org>
To: glibc-cvs@sourceware.org
Subject: [glibc] math: Fix inaccuracy of j0f for x >= 2^127 when sin(x)+cos(x)
 is tiny
X-Act-Checkin: glibc
X-Git-Author: Paul Zimmermann <Paul.Zimmermann@inria.fr>
X-Git-Refname: refs/heads/master
X-Git-Oldrev: 1cfb4715288845ebc55ad664421b48b32de9599c
X-Git-Newrev: b7dd366dbeeff3016d0554c08a5f5c8cfe4c31d1
Message-Id: <20200807193334.26B703857C55@sourceware.org>
Date: Fri,  7 Aug 2020 19:33:34 +0000 (GMT)
X-BeenThere: glibc-cvs@sourceware.org
X-Mailman-Version: 2.1.29
Precedence: list
List-Id: Glibc-cvs mailing list <glibc-cvs.sourceware.org>
List-Unsubscribe: <https://sourceware.org/mailman/options/glibc-cvs>,
 <mailto:glibc-cvs-request@sourceware.org?subject=unsubscribe>
List-Archive: <https://sourceware.org/pipermail/glibc-cvs/>
List-Help: <mailto:glibc-cvs-request@sourceware.org?subject=help>
List-Subscribe: <https://sourceware.org/mailman/listinfo/glibc-cvs>,
 <mailto:glibc-cvs-request@sourceware.org?subject=subscribe>
X-List-Received-Date: Fri, 07 Aug 2020 19:33:34 -0000

https://sourceware.org/git/gitweb.cgi?p=glibc.git;h=b7dd366dbeeff3016d0554c08a5f5c8cfe4c31d1

commit b7dd366dbeeff3016d0554c08a5f5c8cfe4c31d1
Author: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Date:   Fri Aug 7 16:14:53 2020 -0300

    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.

Diff:
---
 math/auto-libm-test-in         |  2 ++
 sysdeps/ieee754/flt-32/e_j0f.c | 17 ++++++++++++++++-
 2 files changed, 18 insertions(+), 1 deletion(-)

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..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)