[PATCH] Fixes for asin[fl].

Steve Kargl sgk at troutmask.apl.washington.edu
Mon Apr 17 00:18:23 UTC 2017


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 (expt<BIAS-1) {	/* |x|<0.5 */
 	    if(expt<ASIN_LINEAR) {	/* if |x| is small, asinl(x)=x */
-		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 = 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


More information about the freebsd-hackers mailing list