[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-numerics
mailing list