[PATCH] Fixes for asin[fl].

[Steve Kargl]

On Mon, Apr 17, 2017 at 02:07:57PM +1000, Bruce Evans wrote:
> 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).

fdlibm will never be updated.  So whatever fdlibm does/did 
is irrelevant.  For e_asinf.c, like most cygnus attributed
code, the conversion from double to float is, ahem, rather
questionable.  As for e_asinl.c, I have doubts to its 
suitability to ld128.

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

Most compilers for the last 20 years know how to convert 1 to
float, double, and long double.  This is in line with getting
rid of (float)0.5 type code. 


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

I'll get to that problem later.  See my sinpi[fl] inplementations.
One does not need a rational approximation (at least with sinpi).

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

If you look at my patch, you'll see that the sqrt() in the 
float version has been replaced with sqrtf().  The float version
cheats by using double in a few places.  I think I know how to
fix that issue.

I also think that sqrt[fl] can be replaced by an in-line specialized
sqrt to avoid the overhead of sqrt[fl].

-- 
Steve
20161221 https://www.youtube.com/watch?v=IbCHE-hONow