From owner-freebsd-hackers@freebsd.org Mon Apr 17 04:08:11 2017 Return-Path: Delivered-To: freebsd-hackers@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 F1CC8D41623; Mon, 17 Apr 2017 04:08:11 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail110.syd.optusnet.com.au (mail110.syd.optusnet.com.au [211.29.132.97]) by mx1.freebsd.org (Postfix) with ESMTP id 767A487B; Mon, 17 Apr 2017 04:08:10 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from besplex.bde.org (c122-106-153-191.carlnfd1.nsw.optusnet.com.au [122.106.153.191]) by mail110.syd.optusnet.com.au (Postfix) with ESMTPS id C1FE010507A; Mon, 17 Apr 2017 14:08:02 +1000 (AEST) Date: Mon, 17 Apr 2017 14:07:57 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Steve Kargl cc: freebsd-numerics@freebsd.org, freebsd-hackers@freebsd.org, freebsd-standards@freebsd.org Subject: Re: [PATCH] Fixes for asin[fl]. In-Reply-To: <20170417001816.GA10277@troutmask.apl.washington.edu> Message-ID: <20170417134317.E931@besplex.bde.org> References: <20170417001816.GA10277@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.2 cv=AYLBJzfG c=1 sm=1 tr=0 a=Tj3pCpwHnMupdyZSltBt7Q==:117 a=Tj3pCpwHnMupdyZSltBt7Q==:17 a=kj9zAlcOel0A:10 a=bSSs6WrcAUF_AD4-SooA:9 a=CjuIK1q_8ugA:10 X-Mailman-Approved-At: Mon, 17 Apr 2017 04:40:20 +0000 X-BeenThere: freebsd-hackers@freebsd.org X-Mailman-Version: 2.1.23 Precedence: list List-Id: Technical Discussions relating to FreeBSD List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 17 Apr 2017 04:08:12 -0000 On Sun, 16 Apr 2017, Steve Kargl wrote: > 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. This is mostly backwards. Don't make style changes away from fdlibm. fdlibm only has e_asin.c, and cygnus' e_asinf.c might have gratuitous differences, and our e_asinl.c is too different due to its hiding of the polynomial coeeficients and evaluation in invtrig.[ch] (there are many style bugs and small errors in the coefficients hidden in the latter). > 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 */ Spelling 1 as "one" is standard style in fdlibm. If fdlibm supported float precision, then it whould do this even more, since this allows spelling all the float constants without an F, as needed to avoid them being double constants. > @@ -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; Removing t is correct. > @@ -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; It was bad to spell 1 as one, but 2 as 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))); This still uses a slow rational function approximation and evaluates it slowly using Horner's method. I have only fixed this in my float version. The polynomial has terms out to S11, and would need too many terms for long doubles. But so does the rational function. > 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; The rational function is not all that slow, since there is a sqrt() as well as a division in some cases, but the division limits accuracy. So does the sqrt(). The float version cheats by using double for sqrt(). Bruce