From owner-freebsd-numerics@freebsd.org Mon Apr 17 00:18:23 2017 Return-Path: Delivered-To: freebsd-numerics@mailman.ysv.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) by mailman.ysv.freebsd.org (Postfix) with ESMTP id B0AEAD41E55; Mon, 17 Apr 2017 00:18:23 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client CN "troutmask", Issuer "troutmask" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id 9D64B1EFA; Mon, 17 Apr 2017 00:18:23 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (localhost [127.0.0.1]) by troutmask.apl.washington.edu (8.15.2/8.15.2) with ESMTPS id v3H0IG66010285 (version=TLSv1.2 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO); Sun, 16 Apr 2017 17:18:16 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.15.2/8.15.2/Submit) id v3H0IGwt010284; Sun, 16 Apr 2017 17:18:16 -0700 (PDT) (envelope-from sgk) Date: Sun, 16 Apr 2017 17:18:16 -0700 From: Steve Kargl To: freebsd-numerics@freebsd.org, freebsd-hackers@freebsd.org, freebsd-standards@freebsd.org Subject: [PATCH] Fixes for asin[fl]. Message-ID: <20170417001816.GA10277@troutmask.apl.washington.edu> Reply-To: sgk@troutmask.apl.washington.edu MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline User-Agent: Mutt/1.7.2 (2016-11-26) X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.23 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 17 Apr 2017 00:18:23 -0000 Please commit. * libm/msun/src/e_asin.c: . Use integer literal constants where possible. This reduces diffs between e_asin[fl].c. . Remove a spurious tab in front of a single line comment. . Sort declaration. . Remove initialization of 't' in declaration statement. It's not needed. * libm/msun/src/e_asinf.c: . Use integer literal constants where possible. This reduces diffs between e_asin[fl].c. . Remove a spurious tab in front of a single line comment. . Sort declaration. . Use sqrtf(x) instead of sqrt(x). The additional precision in the result of sqrt(x) is not required. * libm/msun/src/e_asinl.c: . Remove trailing space in Sunsoft copyright statement. . Use integer literal constants where possible. This reduces diffs between e_asin[fl].c. . Remove a spurious tab in front of a single line comment. . Sort declaration. . Remove initialization of 't' in declaration statement. It's not needed. Reviewed by: md5 Index: src/e_asin.c =================================================================== --- src/e_asin.c (revision 1904) +++ src/e_asin.c (working copy) @@ -50,12 +50,13 @@ __FBSDID("$FreeBSD: head/lib/msun/src/e_ #include "math_private.h" static const double -one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ huge = 1.000e+300, pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ -pio4_hi = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */ - /* coefficient for R(x^2) */ +pio4_hi = 7.85398163397448278999e-01; /* 0x3FE921FB, 0x54442D18 */ + +/* coefficient for R(x^2) */ +static const double pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */ pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */ pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */ @@ -70,7 +71,7 @@ qS4 = 7.70381505559019352791e-02; /* 0x double __ieee754_asin(double x) { - double t=0.0,w,p,q,c,r,s; + double c, p, q, r, s, t, w; int32_t hx,ix; GET_HIGH_WORD(hx,x); ix = hx&0x7fffffff; @@ -83,30 +84,30 @@ __ieee754_asin(double x) return (x-x)/(x-x); /* asin(|x|>1) is NaN */ } else if (ix<0x3fe00000) { /* |x|<0.5 */ if(ix<0x3e500000) { /* if |x| < 2**-26 */ - if(huge+x>one) return x;/* return x with inexact if x!=0*/ + if(huge+x>1) return x;/* return x with inexact if x!=0*/ } t = x*x; p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); - q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); + q = 1+t*(qS1+t*(qS2+t*(qS3+t*qS4))); w = p/q; return x+x*w; } /* 1> |x|>= 0.5 */ - w = one-fabs(x); - t = w*0.5; + w = 1-fabs(x); + t = w/2; p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); - q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); + q = 1+t*(qS1+t*(qS2+t*(qS3+t*qS4))); s = sqrt(t); if(ix>=0x3FEF3333) { /* if |x| > 0.975 */ w = p/q; - t = pio2_hi-(2.0*(s+s*w)-pio2_lo); + t = pio2_hi-(2*(s+s*w)-pio2_lo); } else { w = s; SET_LOW_WORD(w,0); c = (t-w*w)/(s+w); r = p/q; - p = 2.0*s*r-(pio2_lo-2.0*c); - q = pio4_hi-2.0*w; + p = 2*s*r-(pio2_lo-2*c); + q = pio4_hi-2*w; t = pio4_hi-(p-q); } if(hx>0) return t; else return -t; Index: src/e_asinf.c =================================================================== --- src/e_asinf.c (revision 1904) +++ src/e_asinf.c (working copy) @@ -20,9 +20,10 @@ __FBSDID("$FreeBSD: head/lib/msun/src/e_ #include "math_private.h" static const float -one = 1.0000000000e+00, /* 0x3F800000 */ -huge = 1.000e+30, - /* coefficient for R(x^2) */ +huge = 1.000e+30; + +/* coefficient for R(x^2) */ +static const float pS0 = 1.6666586697e-01, pS1 = -4.2743422091e-02, pS2 = -8.6563630030e-03, @@ -35,7 +36,7 @@ float __ieee754_asinf(float x) { double s; - float t,w,p,q; + float p, q, t, w; int32_t hx,ix; GET_FLOAT_WORD(hx,x); ix = hx&0x7fffffff; @@ -45,21 +46,21 @@ __ieee754_asinf(float x) return (x-x)/(x-x); /* asin(|x|>1) is NaN */ } else if (ix<0x3f000000) { /* |x|<0.5 */ if(ix<0x39800000) { /* |x| < 2**-12 */ - if(huge+x>one) return x;/* return x with inexact if x!=0*/ + if(huge+x>1) return x;/* return x with inexact if x!=0*/ } t = x*x; p = t*(pS0+t*(pS1+t*pS2)); - q = one+t*qS1; + q = 1+t*qS1; w = p/q; return x+x*w; } /* 1> |x|>= 0.5 */ - w = one-fabsf(x); - t = w*(float)0.5; + w = 1-fabsf(x); + t = w/2; p = t*(pS0+t*(pS1+t*pS2)); - q = one+t*qS1; - s = sqrt(t); + q = 1+t*qS1; + s = sqrtf(t); w = p/q; - t = pio2-2.0*(s+s*w); + t = pio2-2*(s+s*w); if(hx>0) return t; else return -t; } Index: src/e_asinl.c =================================================================== --- src/e_asinl.c (revision 1904) +++ src/e_asinl.c (working copy) @@ -1,4 +1,3 @@ - /* @(#)e_asin.c 1.3 95/01/18 */ /* FreeBSD: head/lib/msun/src/e_asin.c 176451 2008-02-22 02:30:36Z das */ /* @@ -7,7 +6,7 @@ * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -27,14 +26,13 @@ __FBSDID("$FreeBSD: head/lib/msun/src/e_ #include "math_private.h" static const long double -one = 1.00000000000000000000e+00, huge = 1.000e+300; long double asinl(long double x) { union IEEEl2bits u; - long double t=0.0,w,p,q,c,r,s; + long double c, p, q, r, s, t, w; int16_t expsign, expt; u.e = x; expsign = u.xbits.expsign; @@ -46,7 +44,7 @@ asinl(long double x) return (x-x)/(x-x); /* asin(|x|>1) is NaN */ } else if (exptone) return x;/* return x with inexact if x!=0*/ + if(huge+x>1) return x;/* return x with inexact if x!=0*/ } t = x*x; p = P(t); @@ -55,22 +53,22 @@ asinl(long double x) return x+x*w; } /* 1> |x|>= 0.5 */ - w = one-fabsl(x); - t = w*0.5; + w = 1-fabsl(x); + t = w/2; p = P(t); q = Q(t); s = sqrtl(t); if(u.bits.manh>=THRESH) { /* if |x| is close to 1 */ w = p/q; - t = pio2_hi-(2.0*(s+s*w)-pio2_lo); + t = pio2_hi-(2*(s+s*w)-pio2_lo); } else { u.e = s; u.bits.manl = 0; w = u.e; c = (t-w*w)/(s+w); r = p/q; - p = 2.0*s*r-(pio2_lo-2.0*c); - q = pio4_hi-2.0*w; + p = 2*s*r-(pio2_lo-2*c); + q = pio4_hi-2*w; t = pio4_hi-(p-q); } if(expsign>0) return t; else return -t; -- Steve 20161221 https://www.youtube.com/watch?v=IbCHE-hONow