[cmucl-commit] [git] CMU Common Lisp branch master updated. snapshot-2014-06-52-g7adafd9
Raymond Toy
rtoy at common-lisp.net
Sat Aug 2 06:26:02 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 7adafd921406485dfea0fb8e9290f5ae7f8aa5e5 (commit)
via 1d1ffdf93cb3a67a495b9fe4ea1e3dc679fd401c (commit)
via 21f9b46373c76e44a72b1f7f73cd292397388962 (commit)
via 92c7c5d0c4e9904f1a86a6e3b306ca869d710593 (commit)
from 99afcf7a7ef0b0451cfcb477f8ad241aad930086 (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 7adafd921406485dfea0fb8e9290f5ae7f8aa5e5
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Fri Aug 1 23:25:55 2014 -0700
Use the fdlibm routines for exp and log.
diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index 8f29490..df0bb8c 100644
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -79,8 +79,8 @@
(def-math-rtn "atanh" 1)
;;; Exponential and Logarithmic.
-(def-math-rtn "exp" 1)
-(def-math-rtn "log" 1)
+(def-math-rtn ("__ieee754_exp" %exp) 1)
+(def-math-rtn ("__ieee754_log" %log) 1)
(def-math-rtn "log10" 1)
(def-math-rtn ("__ieee754_pow" %pow) 2)
commit 1d1ffdf93cb3a67a495b9fe4ea1e3dc679fd401c
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Fri Aug 1 23:25:26 2014 -0700
Compile fdlibm routines e_exp.c and e_log.c
diff --git a/src/lisp/Config.x86_darwin b/src/lisp/Config.x86_darwin
index ccf021d..8c7c37b 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
+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
k_sin.o : k_sin.c
$(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<
@@ -44,3 +44,7 @@ s_exmp1.o : s_expm1.c
e_pow.o : e_pow.c
$(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<
+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) $<
commit 21f9b46373c76e44a72b1f7f73cd292397388962
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Fri Aug 1 23:10:46 2014 -0700
Update to use unions.
diff --git a/src/lisp/e_exp.c b/src/lisp/e_exp.c
index e201205..4d94a1e 100644
--- a/src/lisp/e_exp.c
+++ b/src/lisp/e_exp.c
@@ -108,15 +108,17 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
double y,hi,lo,c,t;
int k,xsb;
unsigned hx;
+ union { int i[2]; double d; } ux;
- hx = __HI(x); /* high word of x */
+ ux.d = x;
+ hx = ux.i[HIWORD]; /* high word of x */
xsb = (hx>>31)&1; /* sign bit of x */
hx &= 0x7fffffff; /* high word of |x| */
/* filter out non-finite argument */
if(hx >= 0x40862E42) { /* if |x|>=709.78... */
if(hx>=0x7ff00000) {
- if(((hx&0xfffff)|__LO(x))!=0)
+ if(((hx&0xfffff)|ux.i[LOWORD])!=0)
return x+x; /* NaN */
else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
}
@@ -147,10 +149,14 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
if(k==0) return one-((x*c)/(c-2.0)-x);
else y = one-((lo-(x*c)/(2.0-c))-hi);
if(k >= -1021) {
- __HI(y) += (k<<20); /* add k to y's exponent */
+ ux.d = y;
+ ux.i[HIWORD] += (k<<20); /* add k to y's exponent */
+ y = ux.d;
return y;
} else {
- __HI(y) += ((k+1000)<<20);/* add k to y's exponent */
+ ux.d = y;
+ ux.i[HIWORD] += ((k+1000)<<20);/* add k to y's exponent */
+ y = ux.d;
return y*twom1000;
}
}
diff --git a/src/lisp/e_log.c b/src/lisp/e_log.c
index 3798bc8..4404ce1 100644
--- a/src/lisp/e_log.c
+++ b/src/lisp/e_log.c
@@ -92,9 +92,11 @@ static double zero = 0.0;
double hfsq,f,s,z,R,w,t1,t2,dk;
int k,hx,i,j;
unsigned lx;
+ union { int i[2]; double d; } ux;
- hx = __HI(x); /* high word of x */
- lx = __LO(x); /* low word of x */
+ ux.d = x;
+ hx = ux.i[HIWORD]; /* high word of x */
+ lx = ux.i[LOWORD]; /* low word of x */
k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
@@ -102,13 +104,16 @@ static double zero = 0.0;
return -two54/zero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */
- hx = __HI(x); /* high word of x */
+ ux.d = x;
+ hx = ux.i[HIWORD]; /* high word of x */
}
if (hx >= 0x7ff00000) return x+x;
k += (hx>>20)-1023;
hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000;
- __HI(x) = hx|(i^0x3ff00000); /* normalize x or x/2 */
+ ux.d = x;
+ ux.i[HIWORD] = hx|(i^0x3ff00000); /* normalize x or x/2 */
+ x = ux.d;
k += (i>>20);
f = x-1.0;
if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
commit 92c7c5d0c4e9904f1a86a6e3b306ca869d710593
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Fri Aug 1 23:05:57 2014 -0700
Add fdlibm routines e_exp and e_log, as is.
diff --git a/src/lisp/e_exp.c b/src/lisp/e_exp.c
new file mode 100644
index 0000000..e201205
--- /dev/null
+++ b/src/lisp/e_exp.c
@@ -0,0 +1,156 @@
+
+/* @(#)e_exp.c 1.6 04/04/22 */
+/*
+ * ====================================================
+ * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/* __ieee754_exp(x)
+ * Returns the exponential of x.
+ *
+ * Method
+ * 1. Argument reduction:
+ * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
+ * Given x, find r and integer k such that
+ *
+ * x = k*ln2 + r, |r| <= 0.5*ln2.
+ *
+ * Here r will be represented as r = hi-lo for better
+ * accuracy.
+ *
+ * 2. Approximation of exp(r) by a special rational function on
+ * the interval [0,0.34658]:
+ * Write
+ * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
+ * We use a special Remes algorithm on [0,0.34658] to generate
+ * a polynomial of degree 5 to approximate R. The maximum error
+ * of this polynomial approximation is bounded by 2**-59. In
+ * other words,
+ * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
+ * (where z=r*r, and the values of P1 to P5 are listed below)
+ * and
+ * | 5 | -59
+ * | 2.0+P1*z+...+P5*z - R(z) | <= 2
+ * | |
+ * The computation of exp(r) thus becomes
+ * 2*r
+ * exp(r) = 1 + -------
+ * R - r
+ * r*R1(r)
+ * = 1 + r + ----------- (for better accuracy)
+ * 2 - R1(r)
+ * where
+ * 2 4 10
+ * R1(r) = r - (P1*r + P2*r + ... + P5*r ).
+ *
+ * 3. Scale back to obtain exp(x):
+ * From step 1, we have
+ * exp(x) = 2^k * exp(r)
+ *
+ * Special cases:
+ * exp(INF) is INF, exp(NaN) is NaN;
+ * exp(-INF) is 0, and
+ * for finite argument, only exp(0)=1 is exact.
+ *
+ * Accuracy:
+ * according to an error analysis, the error is always less than
+ * 1 ulp (unit in the last place).
+ *
+ * Misc. info.
+ * For IEEE double
+ * if x > 7.09782712893383973096e+02 then exp(x) overflow
+ * if x < -7.45133219101941108420e+02 then exp(x) underflow
+ *
+ * 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
+#else
+static double
+#endif
+one = 1.0,
+halF[2] = {0.5,-0.5,},
+huge = 1.0e+300,
+twom1000= 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/
+o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
+u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */
+ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
+ -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */
+ln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
+ -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */
+invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
+P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
+P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
+P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
+P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
+P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
+
+
+#ifdef __STDC__
+ double __ieee754_exp(double x) /* default IEEE double exp */
+#else
+ double __ieee754_exp(x) /* default IEEE double exp */
+ double x;
+#endif
+{
+ double y,hi,lo,c,t;
+ int k,xsb;
+ unsigned hx;
+
+ hx = __HI(x); /* high word of x */
+ xsb = (hx>>31)&1; /* sign bit of x */
+ hx &= 0x7fffffff; /* high word of |x| */
+
+ /* filter out non-finite argument */
+ if(hx >= 0x40862E42) { /* if |x|>=709.78... */
+ if(hx>=0x7ff00000) {
+ if(((hx&0xfffff)|__LO(x))!=0)
+ return x+x; /* NaN */
+ else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
+ }
+ if(x > o_threshold) return huge*huge; /* overflow */
+ if(x < u_threshold) return twom1000*twom1000; /* underflow */
+ }
+
+ /* argument reduction */
+ if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
+ if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
+ hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
+ } else {
+ k = (int)(invln2*x+halF[xsb]);
+ t = k;
+ hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */
+ lo = t*ln2LO[0];
+ }
+ x = hi - lo;
+ }
+ else if(hx < 0x3e300000) { /* when |x|<2**-28 */
+ if(huge+x>one) return one+x;/* trigger inexact */
+ }
+ else k = 0;
+
+ /* x is now in primary range */
+ t = x*x;
+ c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
+ if(k==0) return one-((x*c)/(c-2.0)-x);
+ else y = one-((lo-(x*c)/(2.0-c))-hi);
+ if(k >= -1021) {
+ __HI(y) += (k<<20); /* add k to y's exponent */
+ return y;
+ } else {
+ __HI(y) += ((k+1000)<<20);/* add k to y's exponent */
+ return y*twom1000;
+ }
+}
diff --git a/src/lisp/e_log.c b/src/lisp/e_log.c
new file mode 100644
index 0000000..3798bc8
--- /dev/null
+++ b/src/lisp/e_log.c
@@ -0,0 +1,139 @@
+
+/* @(#)e_log.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_log(x)
+ * Return the logrithm of x
+ *
+ * Method :
+ * 1. Argument Reduction: find k and f such that
+ * x = 2^k * (1+f),
+ * where sqrt(2)/2 < 1+f < sqrt(2) .
+ *
+ * 2. Approximation of log(1+f).
+ * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+ * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+ * = 2s + s*R
+ * We use a special Reme algorithm on [0,0.1716] to generate
+ * a polynomial of degree 14 to approximate R The maximum error
+ * of this polynomial approximation is bounded by 2**-58.45. In
+ * other words,
+ * 2 4 6 8 10 12 14
+ * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
+ * (the values of Lg1 to Lg7 are listed in the program)
+ * and
+ * | 2 14 | -58.45
+ * | Lg1*s +...+Lg7*s - R(z) | <= 2
+ * | |
+ * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+ * In order to guarantee error in log below 1ulp, we compute log
+ * by
+ * log(1+f) = f - s*(f - R) (if f is not too large)
+ * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
+ *
+ * 3. Finally, log(x) = k*ln2 + log(1+f).
+ * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+ * Here ln2 is split into two floating point number:
+ * ln2_hi + ln2_lo,
+ * where n*ln2_hi is always exact for |n| < 2000.
+ *
+ * Special cases:
+ * log(x) is NaN with signal if x < 0 (including -INF) ;
+ * log(+INF) is +INF; log(0) is -INF with signal;
+ * log(NaN) is that NaN with no signal.
+ *
+ * Accuracy:
+ * according to an error analysis, the error is always less than
+ * 1 ulp (unit in the last place).
+ *
+ * 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
+#else
+static double
+#endif
+ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
+ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
+two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
+Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
+Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
+Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
+Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
+Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
+Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
+Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
+
+static double zero = 0.0;
+
+#ifdef __STDC__
+ double __ieee754_log(double x)
+#else
+ double __ieee754_log(x)
+ double x;
+#endif
+{
+ double hfsq,f,s,z,R,w,t1,t2,dk;
+ int k,hx,i,j;
+ unsigned lx;
+
+ hx = __HI(x); /* high word of x */
+ lx = __LO(x); /* low word of x */
+
+ k=0;
+ if (hx < 0x00100000) { /* x < 2**-1022 */
+ if (((hx&0x7fffffff)|lx)==0)
+ return -two54/zero; /* log(+-0)=-inf */
+ if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
+ k -= 54; x *= two54; /* subnormal number, scale up x */
+ hx = __HI(x); /* high word of x */
+ }
+ if (hx >= 0x7ff00000) return x+x;
+ k += (hx>>20)-1023;
+ hx &= 0x000fffff;
+ i = (hx+0x95f64)&0x100000;
+ __HI(x) = hx|(i^0x3ff00000); /* normalize x or x/2 */
+ k += (i>>20);
+ f = x-1.0;
+ if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
+ if(f==zero) if(k==0) return zero; else {dk=(double)k;
+ return dk*ln2_hi+dk*ln2_lo;}
+ R = f*f*(0.5-0.33333333333333333*f);
+ if(k==0) return f-R; else {dk=(double)k;
+ return dk*ln2_hi-((R-dk*ln2_lo)-f);}
+ }
+ s = f/(2.0+f);
+ dk = (double)k;
+ z = s*s;
+ i = hx-0x6147a;
+ w = z*z;
+ j = 0x6b851-hx;
+ t1= w*(Lg2+w*(Lg4+w*Lg6));
+ t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
+ i |= j;
+ R = t2+t1;
+ if(i>0) {
+ hfsq=0.5*f*f;
+ if(k==0) return f-(hfsq-s*(hfsq+R)); else
+ return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
+ } else {
+ if(k==0) return f-s*(f-R); else
+ return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
+ }
+}
-----------------------------------------------------------------------
Summary of changes:
src/code/irrat.lisp | 4 +-
src/lisp/Config.x86_darwin | 6 +-
src/lisp/e_exp.c | 162 ++++++++++++++++++++++++++++++++++++++++++++
src/lisp/e_log.c | 144 +++++++++++++++++++++++++++++++++++++++
4 files changed, 313 insertions(+), 3 deletions(-)
create mode 100644 src/lisp/e_exp.c
create mode 100644 src/lisp/e_log.c
hooks/post-receive
--
CMU Common Lisp
More information about the cmucl-commit
mailing list