[cmucl-commit] [git] CMU Common Lisp branch master updated. snapshot-2014-06-56-gc2e152b

Raymond Toy rtoy at common-lisp.net
Sat Aug 2 20:31:36 UTC 2014


This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "CMU Common Lisp".

The branch, master has been updated
       via  c2e152be334f9e1272db737e2eb5056e67def8d0 (commit)
       via  381ee3ea819551345913ff2fe761571d4100aff1 (commit)
       via  9adcb02d4500ccf87a18f8c3b9a9ebbddb4751d0 (commit)
       via  359d84fb9e0a22ae9d6def2bc1a002542354c50c (commit)
      from  7adafd921406485dfea0fb8e9290f5ae7f8aa5e5 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit c2e152be334f9e1272db737e2eb5056e67def8d0
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Sat Aug 2 13:31:18 2014 -0700

    Use the fdlibm asin, acos, and atan routines.

diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index df0bb8c..ef1628b 100644
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -67,10 +67,10 @@
 (def-math-rtn ("fdlibm_sin" %sin) 1)
 (def-math-rtn ("fdlibm_cos" %cos) 1)
 (def-math-rtn ("fdlibm_tan" %tan) 1)
-(def-math-rtn "atan" 1)
+(def-math-rtn ("fdlibm_atan" %atan) 1)
 (def-math-rtn "atan2" 2)
-(def-math-rtn "asin" 1)
-(def-math-rtn "acos" 1)
+(def-math-rtn ("__ieee754_asin" %asin) 1)
+(def-math-rtn ("__ieee754_acos" %acos) 1)
 (def-math-rtn "sinh" 1)
 (def-math-rtn "cosh" 1)
 (def-math-rtn "tanh" 1)

commit 381ee3ea819551345913ff2fe761571d4100aff1
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Sat Aug 2 13:30:49 2014 -0700

    Compile the asin, acos, and atan routines.

diff --git a/src/lisp/Config.x86_darwin b/src/lisp/Config.x86_darwin
index 8c7c37b..9a2c467 100644
--- a/src/lisp/Config.x86_darwin
+++ b/src/lisp/Config.x86_darwin
@@ -18,7 +18,7 @@ OS_LIBS =
 
 EXEC_FINAL_OBJ = exec-final.o
 
-OS_SRC += k_sin.c k_cos.c k_tan.c s_sin.c s_cos.c s_tan.c sincos.c s_log1p.c s_expm1.c e_pow.c e_exp.c e_log.c
+OS_SRC += k_sin.c k_cos.c k_tan.c s_sin.c s_cos.c s_tan.c sincos.c s_log1p.c s_expm1.c e_pow.c e_exp.c e_log.c e_acos.c e_asin.c s_atan.c
 
 k_sin.o : k_sin.c
 	$(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<
@@ -48,3 +48,10 @@ e_exp.o : e_exp.c
 	$(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<
 e_log.o : e_log.c
 	$(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<
+
+e_acos.o : e_acos.c
+	$(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<
+e_asin.o : e_asin.c
+	$(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<
+s_atan.o : s_atan.c
+	$(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<

commit 9adcb02d4500ccf87a18f8c3b9a9ebbddb4751d0
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Sat Aug 2 13:30:17 2014 -0700

    Use unions to access the high and low parts of a double.

diff --git a/src/lisp/e_acos.c b/src/lisp/e_acos.c
index d7c9ed2..bd54a04 100644
--- a/src/lisp/e_acos.c
+++ b/src/lisp/e_acos.c
@@ -66,10 +66,13 @@ qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
 {
 	double z,p,q,r,w,s,c,df;
 	int hx,ix;
-	hx = __HI(x);
+	union { int i[2]; double d; } ux;
+
+        ux.d = x;
+	hx = ux.i[HIWORD];
 	ix = hx&0x7fffffff;
 	if(ix>=0x3ff00000) {	/* |x| >= 1 */
-	    if(((ix-0x3ff00000)|__LO(x))==0) {	/* |x|==1 */
+	    if(((ix-0x3ff00000)|ux.i[LOWORD])==0) {	/* |x|==1 */
 		if(hx>0) return 0.0;		/* acos(1) = 0  */
 		else return pi+2.0*pio2_lo;	/* acos(-1)= pi */
 	    }
@@ -94,7 +97,9 @@ qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
 	    z = (one-x)*0.5;
 	    s = sqrt(z);
 	    df = s;
-	    __LO(df) = 0;
+            ux.d = df;
+	    ux.i[LOWORD] = 0;
+            df = ux.d;
 	    c  = (z-df*df)/(s+df);
 	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
 	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
diff --git a/src/lisp/e_asin.c b/src/lisp/e_asin.c
index 8e37e22..9b476a4 100644
--- a/src/lisp/e_asin.c
+++ b/src/lisp/e_asin.c
@@ -75,10 +75,13 @@ qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
 {
 	double t,w,p,q,c,r,s;
 	int hx,ix;
-	hx = __HI(x);
+	union { int i[2]; double d; } ux;
+
+        ux.d = x;
+	hx = ux.i[HIWORD];
 	ix = hx&0x7fffffff;
 	if(ix>= 0x3ff00000) {		/* |x|>= 1 */
-	    if(((ix-0x3ff00000)|__LO(x))==0)
+	    if(((ix-0x3ff00000)|ux.i[LOWORD])==0)
 		    /* asin(1)=+-pi/2 with inexact */
 		return x*pio2_hi+x*pio2_lo;	
 	    return (x-x)/(x-x);		/* asin(|x|>1) is NaN */   
@@ -103,7 +106,9 @@ qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
 	    t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
 	} else {
 	    w  = s;
-	    __LO(w) = 0;
+            ux.d = w;
+	    ux.i[LOWORD] = 0;
+            w = ux.d;
 	    c  = (t-w*w)/(s+w);
 	    r  = p/q;
 	    p  = 2.0*s*r-(pio2_lo-2.0*c);
diff --git a/src/lisp/s_atan.c b/src/lisp/s_atan.c
index 0093eaf..003f057 100644
--- a/src/lisp/s_atan.c
+++ b/src/lisp/s_atan.c
@@ -83,20 +83,22 @@ one   = 1.0,
 huge   = 1.0e300;
 
 #ifdef __STDC__
-	double atan(double x)
+	double fdlibm_atan(double x)
 #else
-	double atan(x)
+	double fdlibm_atan(x)
 	double x;
 #endif
 {
 	double w,s1,s2,z;
 	int ix,hx,id;
+	union { int i[2]; double d; } ux;
 
-	hx = __HI(x);
+        ux.d = x;
+	hx = ux.i[HIWORD];
 	ix = hx&0x7fffffff;
 	if(ix>=0x44100000) {	/* if |x| >= 2^66 */
 	    if(ix>0x7ff00000||
-		(ix==0x7ff00000&&(__LO(x)!=0)))
+		(ix==0x7ff00000&&(ux.i[LOWORD]!=0)))
 		return x+x;		/* NaN */
 	    if(hx>0) return  atanhi[3]+atanlo[3];
 	    else     return -atanhi[3]-atanlo[3];

commit 359d84fb9e0a22ae9d6def2bc1a002542354c50c
Author: Raymond Toy <toy.raymond at gmail.com>
Date:   Sat Aug 2 08:35:39 2014 -0700

    Import inverse trig functions from fdlibm, as is.

diff --git a/src/lisp/e_acos.c b/src/lisp/e_acos.c
new file mode 100644
index 0000000..d7c9ed2
--- /dev/null
+++ b/src/lisp/e_acos.c
@@ -0,0 +1,105 @@
+
+/* @(#)e_acos.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/* __ieee754_acos(x)
+ * Method :                  
+ *	acos(x)  = pi/2 - asin(x)
+ *	acos(-x) = pi/2 + asin(x)
+ * For |x|<=0.5
+ *	acos(x) = pi/2 - (x + x*x^2*R(x^2))	(see asin.c)
+ * For x>0.5
+ * 	acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
+ *		= 2asin(sqrt((1-x)/2))  
+ *		= 2s + 2s*z*R(z) 	...z=(1-x)/2, s=sqrt(z)
+ *		= 2f + (2c + 2s*z*R(z))
+ *     where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
+ *     for f so that f+c ~ sqrt(z).
+ * For x<-0.5
+ *	acos(x) = pi - 2asin(sqrt((1-|x|)/2))
+ *		= pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z)
+ *
+ * Special cases:
+ *	if x is NaN, return x itself;
+ *	if |x|>1, return NaN with invalid signal.
+ *
+ * Function needed: sqrt
+ */
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+static const double 
+#else
+static double 
+#endif
+one=  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+pi =  3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
+pio2_hi =  1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
+pio2_lo =  6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
+pS0 =  1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
+pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
+pS2 =  2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
+pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
+pS4 =  7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
+pS5 =  3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
+qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
+qS2 =  2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
+qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
+qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
+
+#ifdef __STDC__
+	double __ieee754_acos(double x)
+#else
+	double __ieee754_acos(x)
+	double x;
+#endif
+{
+	double z,p,q,r,w,s,c,df;
+	int hx,ix;
+	hx = __HI(x);
+	ix = hx&0x7fffffff;
+	if(ix>=0x3ff00000) {	/* |x| >= 1 */
+	    if(((ix-0x3ff00000)|__LO(x))==0) {	/* |x|==1 */
+		if(hx>0) return 0.0;		/* acos(1) = 0  */
+		else return pi+2.0*pio2_lo;	/* acos(-1)= pi */
+	    }
+	    return (x-x)/(x-x);		/* acos(|x|>1) is NaN */
+	}
+	if(ix<0x3fe00000) {	/* |x| < 0.5 */
+	    if(ix<=0x3c600000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/
+	    z = x*x;
+	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+	    r = p/q;
+	    return pio2_hi - (x - (pio2_lo-x*r));
+	} else  if (hx<0) {		/* x < -0.5 */
+	    z = (one+x)*0.5;
+	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+	    s = sqrt(z);
+	    r = p/q;
+	    w = r*s-pio2_lo;
+	    return pi - 2.0*(s+w);
+	} else {			/* x > 0.5 */
+	    z = (one-x)*0.5;
+	    s = sqrt(z);
+	    df = s;
+	    __LO(df) = 0;
+	    c  = (z-df*df)/(s+df);
+	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+	    r = p/q;
+	    w = r*s+c;
+	    return 2.0*(df+w);
+	}
+}
diff --git a/src/lisp/e_asin.c b/src/lisp/e_asin.c
new file mode 100644
index 0000000..8e37e22
--- /dev/null
+++ b/src/lisp/e_asin.c
@@ -0,0 +1,114 @@
+
+/* @(#)e_asin.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/* __ieee754_asin(x)
+ * Method :                  
+ *	Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
+ *	we approximate asin(x) on [0,0.5] by
+ *		asin(x) = x + x*x^2*R(x^2)
+ *	where
+ *		R(x^2) is a rational approximation of (asin(x)-x)/x^3 
+ *	and its remez error is bounded by
+ *		|(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
+ *
+ *	For x in [0.5,1]
+ *		asin(x) = pi/2-2*asin(sqrt((1-x)/2))
+ *	Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
+ *	then for x>0.98
+ *		asin(x) = pi/2 - 2*(s+s*z*R(z))
+ *			= pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
+ *	For x<=0.98, let pio4_hi = pio2_hi/2, then
+ *		f = hi part of s;
+ *		c = sqrt(z) - f = (z-f*f)/(s+f) 	...f+c=sqrt(z)
+ *	and
+ *		asin(x) = pi/2 - 2*(s+s*z*R(z))
+ *			= pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
+ *			= pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
+ *
+ * Special cases:
+ *	if x is NaN, return x itself;
+ *	if |x|>1, return NaN with invalid signal.
+ *
+ */
+
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+static const double 
+#else
+static double 
+#endif
+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) */
+pS0 =  1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
+pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
+pS2 =  2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
+pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
+pS4 =  7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
+pS5 =  3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
+qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
+qS2 =  2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
+qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
+qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
+
+#ifdef __STDC__
+	double __ieee754_asin(double x)
+#else
+	double __ieee754_asin(x)
+	double x;
+#endif
+{
+	double t,w,p,q,c,r,s;
+	int hx,ix;
+	hx = __HI(x);
+	ix = hx&0x7fffffff;
+	if(ix>= 0x3ff00000) {		/* |x|>= 1 */
+	    if(((ix-0x3ff00000)|__LO(x))==0)
+		    /* asin(1)=+-pi/2 with inexact */
+		return x*pio2_hi+x*pio2_lo;	
+	    return (x-x)/(x-x);		/* asin(|x|>1) is NaN */   
+	} else if (ix<0x3fe00000) {	/* |x|<0.5 */
+	    if(ix<0x3e400000) {		/* if |x| < 2**-27 */
+		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)));
+		w = p/q;
+		return x+x*w;
+	}
+	/* 1> |x|>= 0.5 */
+	w = one-fabs(x);
+	t = w*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 = sqrt(t);
+	if(ix>=0x3FEF3333) { 	/* if |x| > 0.975 */
+	    w = p/q;
+	    t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
+	} else {
+	    w  = s;
+	    __LO(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;
+	    t  = pio4_hi-(p-q);
+	}    
+	if(hx>0) return t; else return -t;    
+}
diff --git a/src/lisp/s_atan.c b/src/lisp/s_atan.c
new file mode 100644
index 0000000..0093eaf
--- /dev/null
+++ b/src/lisp/s_atan.c
@@ -0,0 +1,134 @@
+
+/* @(#)s_atan.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ */
+
+/* atan(x)
+ * Method
+ *   1. Reduce x to positive by atan(x) = -atan(-x).
+ *   2. According to the integer k=4t+0.25 chopped, t=x, the argument
+ *      is further reduced to one of the following intervals and the
+ *      arctangent of t is evaluated by the corresponding formula:
+ *
+ *      [0,7/16]      atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
+ *      [7/16,11/16]  atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
+ *      [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
+ *      [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
+ *      [39/16,INF]   atan(x) = atan(INF) + atan( -1/t )
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following 
+ * constants. The decimal values may be used, provided that the 
+ * compiler will convert from decimal to binary accurately enough 
+ * to produce the hexadecimal values shown.
+ */
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+static const double atanhi[] = {
+#else
+static double atanhi[] = {
+#endif
+  4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
+  7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
+  9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
+  1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
+};
+
+#ifdef __STDC__
+static const double atanlo[] = {
+#else
+static double atanlo[] = {
+#endif
+  2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
+  3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
+  1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
+  6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
+};
+
+#ifdef __STDC__
+static const double aT[] = {
+#else
+static double aT[] = {
+#endif
+  3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
+ -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
+  1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
+ -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
+  9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */
+ -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
+  6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */
+ -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
+  4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */
+ -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
+  1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
+};
+
+#ifdef __STDC__
+	static const double 
+#else
+	static double 
+#endif
+one   = 1.0,
+huge   = 1.0e300;
+
+#ifdef __STDC__
+	double atan(double x)
+#else
+	double atan(x)
+	double x;
+#endif
+{
+	double w,s1,s2,z;
+	int ix,hx,id;
+
+	hx = __HI(x);
+	ix = hx&0x7fffffff;
+	if(ix>=0x44100000) {	/* if |x| >= 2^66 */
+	    if(ix>0x7ff00000||
+		(ix==0x7ff00000&&(__LO(x)!=0)))
+		return x+x;		/* NaN */
+	    if(hx>0) return  atanhi[3]+atanlo[3];
+	    else     return -atanhi[3]-atanlo[3];
+	} if (ix < 0x3fdc0000) {	/* |x| < 0.4375 */
+	    if (ix < 0x3e200000) {	/* |x| < 2^-29 */
+		if(huge+x>one) return x;	/* raise inexact */
+	    }
+	    id = -1;
+	} else {
+	x = fabs(x);
+	if (ix < 0x3ff30000) {		/* |x| < 1.1875 */
+	    if (ix < 0x3fe60000) {	/* 7/16 <=|x|<11/16 */
+		id = 0; x = (2.0*x-one)/(2.0+x); 
+	    } else {			/* 11/16<=|x|< 19/16 */
+		id = 1; x  = (x-one)/(x+one); 
+	    }
+	} else {
+	    if (ix < 0x40038000) {	/* |x| < 2.4375 */
+		id = 2; x  = (x-1.5)/(one+1.5*x);
+	    } else {			/* 2.4375 <= |x| < 2^66 */
+		id = 3; x  = -1.0/x;
+	    }
+	}}
+    /* end of argument reduction */
+	z = x*x;
+	w = z*z;
+    /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
+	s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
+	s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
+	if (id<0) return x - x*(s1+s2);
+	else {
+	    z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
+	    return (hx<0)? -z:z;
+	}
+}

-----------------------------------------------------------------------

Summary of changes:
 src/code/irrat.lisp        |    6 +-
 src/lisp/Config.x86_darwin |    9 ++-
 src/lisp/e_acos.c          |  110 +++++++++++++++++++++++++++++++++++
 src/lisp/e_asin.c          |  119 ++++++++++++++++++++++++++++++++++++++
 src/lisp/s_atan.c          |  136 ++++++++++++++++++++++++++++++++++++++++++++
 5 files changed, 376 insertions(+), 4 deletions(-)
 create mode 100644 src/lisp/e_acos.c
 create mode 100644 src/lisp/e_asin.c
 create mode 100644 src/lisp/s_atan.c


hooks/post-receive
-- 
CMU Common Lisp


More information about the cmucl-commit mailing list