[cmucl-commit] [git] CMU Common Lisp branch master updated. snapshot-2014-06-40-ge0b1f9f
Raymond Toy
rtoy at common-lisp.net
Sat Aug 2 02:31:53 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 e0b1f9f8b2142397cbf4ea76dd3ba862862baa49 (commit)
via e19642f9bb877846f2fa05701bebcef1073821f3 (commit)
via ae8eb31e7757d351eca8fd5b9693637843cda611 (commit)
from 6795c643cd67e3ea2d950aa249cf39f929dc7798 (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 e0b1f9f8b2142397cbf4ea76dd3ba862862baa49
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Fri Aug 1 19:31:39 2014 -0700
Compile s_log1p and s_expm1.
diff --git a/src/lisp/Config.x86_darwin b/src/lisp/Config.x86_darwin
index 3dcdc61..70c2df4 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
+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
k_sin.o : k_sin.c
$(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<
@@ -37,3 +37,7 @@ s_tan.o : s_tan.c
sincos.o : sincos.c
$(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<
+s_log1p.o : s_log1p.c
+ $(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<
+s_exmp1.o : s_expm1.c
+ $(CC) -c $(CFLAGS) $(CPPFLAGS) $(CC_REM_PIO2) $<
commit e19642f9bb877846f2fa05701bebcef1073821f3
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Fri Aug 1 19:30:45 2014 -0700
Modify to use unions and change name to include fdlibm prefix.
diff --git a/src/lisp/s_expm1.c b/src/lisp/s_expm1.c
index e4dd820..d6a1827 100644
--- a/src/lisp/s_expm1.c
+++ b/src/lisp/s_expm1.c
@@ -127,17 +127,20 @@ Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */
Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
#ifdef __STDC__
- double expm1(double x)
+ double fdlibm_expm1(double x)
#else
- double expm1(x)
+ double fdlibm_expm1(x)
double x;
#endif
{
double y,hi,lo,c,t,e,hxs,hfx,r1;
int k,xsb;
unsigned hx;
+ union { int i[2]; double d; } ux;
+ union { int i[2]; double d; } utmp;
- hx = __HI(x); /* high word of x */
+ ux.d = x;
+ hx = ux.i[HIWORD]; /* high word of x */
xsb = hx&0x80000000; /* sign bit of x */
if(xsb==0) y=x; else y= -x; /* y = |x| */
hx &= 0x7fffffff; /* high word of |x| */
@@ -146,7 +149,7 @@ Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */
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:-1.0;/* exp(+-inf)={inf,-1} */
}
@@ -196,19 +199,29 @@ Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
else return one+2.0*(x-e);
if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */
y = one-(e-x);
- __HI(y) += (k<<20); /* add k to y's exponent */
+ utmp.d = y;
+ utmp.i[HIWORD] += (k<<20); /* add k to y's exponent */
+ y = utmp.d;
return y-one;
}
t = one;
if(k<20) {
- __HI(t) = 0x3ff00000 - (0x200000>>k); /* t=1-2^-k */
+ utmp.d = t;
+ utmp.i[HIWORD] = 0x3ff00000 - (0x200000>>k); /* t=1-2^-k */
+ t = utmp.d;
y = t-(e-x);
- __HI(y) += (k<<20); /* add k to y's exponent */
+ utmp.d = y;
+ utmp.i[HIWORD] += (k<<20); /* add k to y's exponent */
+ y = utmp.d;
} else {
- __HI(t) = ((0x3ff-k)<<20); /* 2^-k */
+ utmp.d = t;
+ utmp.i[HIWORD] = ((0x3ff-k)<<20); /* 2^-k */
+ t = utmp.d;
y = x-(e+t);
y += one;
- __HI(y) += (k<<20); /* add k to y's exponent */
+ utmp.d = y;
+ utmp.i[HIWORD] += (k<<20); /* add k to y's exponent */
+ y = utmp.d;
}
}
return y;
diff --git a/src/lisp/s_log1p.c b/src/lisp/s_log1p.c
index 3d1059a..eaec38a 100644
--- a/src/lisp/s_log1p.c
+++ b/src/lisp/s_log1p.c
@@ -97,16 +97,18 @@ Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static double zero = 0.0;
#ifdef __STDC__
- double log1p(double x)
+ double fdlibm_log1p(double x)
#else
- double log1p(x)
+ double fdlibm_log1p(x)
double x;
#endif
{
double hfsq,f,c,s,z,R,u;
int k,hx,hu,ax;
+ 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 */
ax = hx&0x7fffffff;
k = 1;
@@ -128,23 +130,29 @@ static double zero = 0.0;
if (hx >= 0x7ff00000) return x+x;
if(k!=0) {
if(hx<0x43400000) {
- u = 1.0+x;
- hu = __HI(u); /* high word of u */
+ u = 1.0+x;
+ ux.d = u;
+ hu = ux.i[HIWORD]; /* high word of u */
k = (hu>>20)-1023;
c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */
c /= u;
} else {
u = x;
- hu = __HI(u); /* high word of u */
+ ux.d = u;
+ hu = ux.i[HIWORD]; /* high word of u */
k = (hu>>20)-1023;
c = 0;
}
hu &= 0x000fffff;
if(hu<0x6a09e) {
- __HI(u) = hu|0x3ff00000; /* normalize u */
+ ux.d = u;
+ ux.i[HIWORD] = hu|0x3ff00000; /* normalize u */
+ u = ux.d;
} else {
- k += 1;
- __HI(u) = hu|0x3fe00000; /* normalize u/2 */
+ k += 1;
+ ux.d = u;
+ ux.i[HIWORD] = hu|0x3fe00000; /* normalize u/2 */
+ u = ux.d;
hu = (0x00100000-hu)>>2;
}
f = u-1.0;
commit ae8eb31e7757d351eca8fd5b9693637843cda611
Author: Raymond Toy <toy.raymond at gmail.com>
Date: Fri Aug 1 19:21:05 2014 -0700
Import expm1 from fdlibm, as is.
diff --git a/src/lisp/s_expm1.c b/src/lisp/s_expm1.c
new file mode 100644
index 0000000..e4dd820
--- /dev/null
+++ b/src/lisp/s_expm1.c
@@ -0,0 +1,215 @@
+
+/* @(#)s_expm1.c 1.5 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.
+ * ====================================================
+ */
+
+/* expm1(x)
+ * Returns exp(x)-1, the exponential of x minus 1.
+ *
+ * Method
+ * 1. Argument reduction:
+ * Given x, find r and integer k such that
+ *
+ * x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658
+ *
+ * Here a correction term c will be computed to compensate
+ * the error in r when rounded to a floating-point number.
+ *
+ * 2. Approximating expm1(r) by a special rational function on
+ * the interval [0,0.34658]:
+ * Since
+ * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
+ * we define R1(r*r) by
+ * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r)
+ * That is,
+ * R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
+ * = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
+ * = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
+ * We use a special Remes algorithm on [0,0.347] to generate
+ * a polynomial of degree 5 in r*r to approximate R1. The
+ * maximum error of this polynomial approximation is bounded
+ * by 2**-61. In other words,
+ * R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
+ * where Q1 = -1.6666666666666567384E-2,
+ * Q2 = 3.9682539681370365873E-4,
+ * Q3 = -9.9206344733435987357E-6,
+ * Q4 = 2.5051361420808517002E-7,
+ * Q5 = -6.2843505682382617102E-9;
+ * (where z=r*r, and the values of Q1 to Q5 are listed below)
+ * with error bounded by
+ * | 5 | -61
+ * | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2
+ * | |
+ *
+ * expm1(r) = exp(r)-1 is then computed by the following
+ * specific way which minimize the accumulation rounding error:
+ * 2 3
+ * r r [ 3 - (R1 + R1*r/2) ]
+ * expm1(r) = r + --- + --- * [--------------------]
+ * 2 2 [ 6 - r*(3 - R1*r/2) ]
+ *
+ * To compensate the error in the argument reduction, we use
+ * expm1(r+c) = expm1(r) + c + expm1(r)*c
+ * ~ expm1(r) + c + r*c
+ * Thus c+r*c will be added in as the correction terms for
+ * expm1(r+c). Now rearrange the term to avoid optimization
+ * screw up:
+ * ( 2 2 )
+ * ({ ( r [ R1 - (3 - R1*r/2) ] ) } r )
+ * expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
+ * ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 )
+ * ( )
+ *
+ * = r - E
+ * 3. Scale back to obtain expm1(x):
+ * From step 1, we have
+ * expm1(x) = either 2^k*[expm1(r)+1] - 1
+ * = or 2^k*[expm1(r) + (1-2^-k)]
+ * 4. Implementation notes:
+ * (A). To save one multiplication, we scale the coefficient Qi
+ * to Qi*2^i, and replace z by (x^2)/2.
+ * (B). To achieve maximum accuracy, we compute expm1(x) by
+ * (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
+ * (ii) if k=0, return r-E
+ * (iii) if k=-1, return 0.5*(r-E)-0.5
+ * (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E)
+ * else return 1.0+2.0*(r-E);
+ * (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
+ * (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else
+ * (vii) return 2^k(1-((E+2^-k)-r))
+ *
+ * Special cases:
+ * expm1(INF) is INF, expm1(NaN) is NaN;
+ * expm1(-INF) is -1, and
+ * for finite argument, only expm1(0)=0 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 expm1(x) overflow
+ *
+ * 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,
+huge = 1.0e+300,
+tiny = 1.0e-300,
+o_threshold = 7.09782712893383973096e+02,/* 0x40862E42, 0xFEFA39EF */
+ln2_hi = 6.93147180369123816490e-01,/* 0x3fe62e42, 0xfee00000 */
+ln2_lo = 1.90821492927058770002e-10,/* 0x3dea39ef, 0x35793c76 */
+invln2 = 1.44269504088896338700e+00,/* 0x3ff71547, 0x652b82fe */
+ /* scaled coefficients related to expm1 */
+Q1 = -3.33333333333331316428e-02, /* BFA11111 111110F4 */
+Q2 = 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */
+Q3 = -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */
+Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */
+Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
+
+#ifdef __STDC__
+ double expm1(double x)
+#else
+ double expm1(x)
+ double x;
+#endif
+{
+ double y,hi,lo,c,t,e,hxs,hfx,r1;
+ int k,xsb;
+ unsigned hx;
+
+ hx = __HI(x); /* high word of x */
+ xsb = hx&0x80000000; /* sign bit of x */
+ if(xsb==0) y=x; else y= -x; /* y = |x| */
+ hx &= 0x7fffffff; /* high word of |x| */
+
+ /* filter out huge and non-finite argument */
+ if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */
+ 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:-1.0;/* exp(+-inf)={inf,-1} */
+ }
+ if(x > o_threshold) return huge*huge; /* overflow */
+ }
+ if(xsb!=0) { /* x < -56*ln2, return -1.0 with inexact */
+ if(x+tiny<0.0) /* raise inexact */
+ return tiny-one; /* return -1 */
+ }
+ }
+
+ /* argument reduction */
+ if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
+ if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
+ if(xsb==0)
+ {hi = x - ln2_hi; lo = ln2_lo; k = 1;}
+ else
+ {hi = x + ln2_hi; lo = -ln2_lo; k = -1;}
+ } else {
+ k = invln2*x+((xsb==0)?0.5:-0.5);
+ t = k;
+ hi = x - t*ln2_hi; /* t*ln2_hi is exact here */
+ lo = t*ln2_lo;
+ }
+ x = hi - lo;
+ c = (hi-x)-lo;
+ }
+ else if(hx < 0x3c900000) { /* when |x|<2**-54, return x */
+ t = huge+x; /* return x with inexact flags when x!=0 */
+ return x - (t-(huge+x));
+ }
+ else k = 0;
+
+ /* x is now in primary range */
+ hfx = 0.5*x;
+ hxs = x*hfx;
+ r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
+ t = 3.0-r1*hfx;
+ e = hxs*((r1-t)/(6.0 - x*t));
+ if(k==0) return x - (x*e-hxs); /* c is 0 */
+ else {
+ e = (x*(e-c)-c);
+ e -= hxs;
+ if(k== -1) return 0.5*(x-e)-0.5;
+ if(k==1)
+ if(x < -0.25) return -2.0*(e-(x+0.5));
+ else return one+2.0*(x-e);
+ if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */
+ y = one-(e-x);
+ __HI(y) += (k<<20); /* add k to y's exponent */
+ return y-one;
+ }
+ t = one;
+ if(k<20) {
+ __HI(t) = 0x3ff00000 - (0x200000>>k); /* t=1-2^-k */
+ y = t-(e-x);
+ __HI(y) += (k<<20); /* add k to y's exponent */
+ } else {
+ __HI(t) = ((0x3ff-k)<<20); /* 2^-k */
+ y = x-(e+t);
+ y += one;
+ __HI(y) += (k<<20); /* add k to y's exponent */
+ }
+ }
+ return y;
+}
-----------------------------------------------------------------------
Summary of changes:
src/lisp/Config.x86_darwin | 6 +-
src/lisp/s_expm1.c | 228 ++++++++++++++++++++++++++++++++++++++++++++
src/lisp/s_log1p.c | 26 +++--
3 files changed, 250 insertions(+), 10 deletions(-)
create mode 100644 src/lisp/s_expm1.c
hooks/post-receive
--
CMU Common Lisp
More information about the cmucl-commit
mailing list