[PATCH] Fixes for asin[fl].

[Bruce Evans]

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