From mboxrd@z Thu Jan 1 00:00:00 1970 X-Msuck: nntp://news.gmane.org/gmane.linux.lib.musl.general/11134 Path: news.gmane.org!.POSTED!not-for-mail From: Szabolcs Nagy Newsgroups: gmane.linux.lib.musl.general Subject: [PATCH] fix threshold constants in j0f, y0f, j1f, y1f Date: Wed, 15 Mar 2017 02:55:49 +0100 Message-ID: <20170315015549.GL2082@port70.net> Reply-To: musl@lists.openwall.com NNTP-Posting-Host: blaine.gmane.org Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii X-Trace: blaine.gmane.org 1489542967 5374 195.159.176.226 (15 Mar 2017 01:56:07 GMT) X-Complaints-To: usenet@blaine.gmane.org NNTP-Posting-Date: Wed, 15 Mar 2017 01:56:07 +0000 (UTC) User-Agent: Mutt/1.6.0 (2016-04-01) To: musl@lists.openwall.com Original-X-From: musl-return-11149-gllmg-musl=m.gmane.org@lists.openwall.com Wed Mar 15 02:56:02 2017 Return-path: Envelope-to: gllmg-musl@m.gmane.org Original-Received: from mother.openwall.net ([195.42.179.200]) by blaine.gmane.org with smtp (Exim 4.84_2) (envelope-from ) id 1cnyAR-0000a6-Jc for gllmg-musl@m.gmane.org; Wed, 15 Mar 2017 02:55:59 +0100 Original-Received: (qmail 11827 invoked by uid 550); 15 Mar 2017 01:56:03 -0000 Mailing-List: contact musl-help@lists.openwall.com; run by ezmlm Precedence: bulk List-Post: List-Help: List-Unsubscribe: List-Subscribe: List-ID: Original-Received: (qmail 11783 invoked from network); 15 Mar 2017 01:56:01 -0000 Mail-Followup-To: musl@lists.openwall.com Content-Disposition: inline Xref: news.gmane.org gmane.linux.lib.musl.general:11134 Archived-At: partly following freebsd rev 279491 https://svnweb.freebsd.org/base?view=revision&revision=279491 (musl had some of the fixes before freebsd). the change should not matter much for j0f, y0f, but it improves j1f and y1f in [2.5,~3.75] (that is [0x40200000,~0x40700000]). near roots (e.g. around 3.8317 for j1f) there are still large ulp errors. dropped code that tried to raise inexact. --- src/math/j0f.c | 8 ++++---- src/math/j1f.c | 17 ++++++++--------- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/src/math/j0f.c b/src/math/j0f.c index 45883dc4..fab554a3 100644 --- a/src/math/j0f.c +++ b/src/math/j0f.c @@ -208,8 +208,8 @@ static float pzerof(float x) GET_FLOAT_WORD(ix, x); ix &= 0x7fffffff; if (ix >= 0x41000000){p = pR8; q = pS8;} - else if (ix >= 0x40f71c58){p = pR5; q = pS5;} - else if (ix >= 0x4036db68){p = pR3; q = pS3;} + else if (ix >= 0x409173eb){p = pR5; q = pS5;} + else if (ix >= 0x4036d917){p = pR3; q = pS3;} else /*ix >= 0x40000000*/ {p = pR2; q = pS2;} z = 1.0f/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); @@ -304,8 +304,8 @@ static float qzerof(float x) GET_FLOAT_WORD(ix, x); ix &= 0x7fffffff; if (ix >= 0x41000000){p = qR8; q = qS8;} - else if (ix >= 0x40f71c58){p = qR5; q = qS5;} - else if (ix >= 0x4036db68){p = qR3; q = qS3;} + else if (ix >= 0x409173eb){p = qR5; q = qS5;} + else if (ix >= 0x4036d917){p = qR3; q = qS3;} else /*ix >= 0x40000000*/ {p = qR2; q = qS2;} z = 1.0f/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); diff --git a/src/math/j1f.c b/src/math/j1f.c index 58875af9..3434c53d 100644 --- a/src/math/j1f.c +++ b/src/math/j1f.c @@ -74,14 +74,13 @@ float j1f(float x) return 1/(x*x); if (ix >= 0x40000000) /* |x| >= 2 */ return common(ix, fabsf(x), 0, sign); - if (ix >= 0x32000000) { /* |x| >= 2**-27 */ + if (ix >= 0x39000000) { /* |x| >= 2**-13 */ z = x*x; r = z*(r00+z*(r01+z*(r02+z*r03))); s = 1+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); z = 0.5f + r/s; } else - /* raise inexact if x!=0 */ - z = 0.5f + x; + z = 0.5f; return z*x; } @@ -114,7 +113,7 @@ float y1f(float x) return 1/x; if (ix >= 0x40000000) /* |x| >= 2.0 */ return common(ix,x,1,0); - if (ix < 0x32000000) /* x < 2**-27 */ + if (ix < 0x33000000) /* x < 2**-25 */ return -tpi/x; z = x*x; u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); @@ -205,8 +204,8 @@ static float ponef(float x) GET_FLOAT_WORD(ix, x); ix &= 0x7fffffff; if (ix >= 0x41000000){p = pr8; q = ps8;} - else if (ix >= 0x40f71c58){p = pr5; q = ps5;} - else if (ix >= 0x4036db68){p = pr3; q = ps3;} + else if (ix >= 0x409173eb){p = pr5; q = ps5;} + else if (ix >= 0x4036d917){p = pr3; q = ps3;} else /*ix >= 0x40000000*/ {p = pr2; q = ps2;} z = 1.0f/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); @@ -300,9 +299,9 @@ static float qonef(float x) GET_FLOAT_WORD(ix, x); ix &= 0x7fffffff; - if (ix >= 0x40200000){p = qr8; q = qs8;} - else if (ix >= 0x40f71c58){p = qr5; q = qs5;} - else if (ix >= 0x4036db68){p = qr3; q = qs3;} + if (ix >= 0x41000000){p = qr8; q = qs8;} + else if (ix >= 0x409173eb){p = qr5; q = qs5;} + else if (ix >= 0x4036d917){p = qr3; q = qs3;} else /*ix >= 0x40000000*/ {p = qr2; q = qs2;} z = 1.0f/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); -- 2.11.0