Fix some problems with asinf(), acosf(), atanf(), and atan2f():

- Adjust several constants for float precision. Some thresholds
  that were appropriate for double precision were never changed
  when these routines were converted to float precision. This
  has an impact on performance but not accuracy. (Submitted by bde.)

- Reduce the degrees of the polynomials used. A smaller degree
  suffices for float precision.

- In asinf(), use double arithmetic in part of the calculation to
  avoid a corner case and some complicated arithmetic involving a
  division and some buggy constants. This improves performance and
  accuracy.

Max error (ulps):
         asinf  acosf  atanf
before   0.925  0.782  0.852
after    0.743  0.804  0.852

As bde points out, it's cheaper for asin*() and acos*() to use
polynomials instead of rational functions, but that's a task for
another day.

1 files changed, 18 insertions, 35 deletions

diff --git a/lib/msun/src/e_asinf.c b/lib/msun/src/e_asinf.c
index f830878..a4c3f00 100644
--- a/lib/msun/src/e_asinf.c
+++ b/lib/msun/src/e_asinf.c
@@ -22,62 +22,45 @@ __FBSDID("$FreeBSD$");
 static const float
 one =  1.0000000000e+00, /* 0x3F800000 */
 huge =  1.000e+30,
-pio2_hi =  1.5707962513e+00, /* 0x3fc90fda */
-pio2_lo =  7.5497894159e-08, /* 0x33a22168 */
-pio4_hi =  7.8539812565e-01, /* 0x3f490fda */
 	/* coefficient for R(x^2) */
-pS0 =  1.6666667163e-01, /* 0x3e2aaaab */
-pS1 = -3.2556581497e-01, /* 0xbea6b090 */
-pS2 =  2.0121252537e-01, /* 0x3e4e0aa8 */
-pS3 = -4.0055535734e-02, /* 0xbd241146 */
-pS4 =  7.9153501429e-04, /* 0x3a4f7f04 */
-pS5 =  3.4793309169e-05, /* 0x3811ef08 */
-qS1 = -2.4033949375e+00, /* 0xc019d139 */
-qS2 =  2.0209457874e+00, /* 0x4001572d */
-qS3 = -6.8828397989e-01, /* 0xbf303361 */
-qS4 =  7.7038154006e-02; /* 0x3d9dc62e */
+pS0 =  1.6666586697e-01,
+pS1 = -4.2743422091e-02,
+pS2 = -8.6563630030e-03,
+qS1 = -7.0662963390e-01;
+
+static const double
+pio2 =  1.570796326794896558e+00;
 
 float
 __ieee754_asinf(float x)
 {
-	float t=0.0,w,p,q,c,r,s;
+	double s;
+	float t=0.0,w,p,q,c,r;
 	int32_t hx,ix;
 	GET_FLOAT_WORD(hx,x);
 	ix = hx&0x7fffffff;
 	if(ix==0x3f800000) {
 		/* asin(1)=+-pi/2 with inexact */
-	    return x*pio2_hi+x*pio2_lo;
+	    return x*pio2;
 	} else if(ix> 0x3f800000) {	/* |x|>= 1 */
 	    return (x-x)/(x-x);		/* asin(|x|>1) is NaN */
 	} else if (ix<0x3f000000) {	/* |x|<0.5 */
-	    if(ix<0x32000000) {		/* if |x| < 2**-27 */
+	    if(ix<0x39800000) {		/* if |x| < 2**-12 */
 		if(huge+x>one) return x;/* return x with inexact if x!=0*/
 	    } else
 		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)));
+	        p = t*(pS0+t*(pS1+t*pS2));
+	        q = one+t*qS1;
 		w = p/q;
 		return x+x*w;
 	}
 	/* 1> |x|>= 0.5 */
 	w = one-fabsf(x);
 	t = w*(float)0.5;
-	p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
-	q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
-	s = __ieee754_sqrtf(t);
-	if(ix>=0x3F79999A) { 	/* if |x| > 0.975 */
-	    w = p/q;
-	    t = pio2_hi-((float)2.0*(s+s*w)-pio2_lo);
-	} else {
-	    int32_t iw;
-	    w  = s;
-	    GET_FLOAT_WORD(iw,w);
-	    SET_FLOAT_WORD(w,iw&0xfffff000);
-	    c  = (t-w*w)/(s+w);
-	    r  = p/q;
-	    p  = (float)2.0*s*r-(pio2_lo-(float)2.0*c);
-	    q  = pio4_hi-(float)2.0*w;
-	    t  = pio4_hi-(p-q);
-	}
+	p = t*(pS0+t*(pS1+t*pS2));
+	q = one+t*qS1;
+	s = __ieee754_sqrt(t);
+	w = p/q;
+	t = pio2-2.0*(s+s*w);
 	if(hx>0) return t; else return -t;
 }