X-Git-Url: http://nsz.repo.hu/git/?a=blobdiff_plain;f=src%2Fmath%2Flog.c;h=e61e113d41af91f86d4482a979c84966be66fe33;hb=8041af59881219c32267c3491bee43591d3c3fe6;hp=98051205f80ec607d416a645a1d77d63cb0eebc8;hpb=97721a5508415a2f10eb068e022093811c9ff8be;p=musl

diff --git a/src/math/log.c b/src/math/log.c
index 98051205..e61e113d 100644
--- a/src/math/log.c
+++ b/src/math/log.c
@@ -10,7 +10,7 @@
  * ====================================================
  */
 /* log(x)
- * Return the logrithm of x
+ * Return the logarithm of x
  *
  * Method :
  *   1. Argument Reduction: find k and f such that
@@ -60,12 +60,12 @@
  * to produce the hexadecimal values shown.
  */
 
-#include "libm.h"
+#include <math.h>
+#include <stdint.h>
 
 static const double
 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 */
@@ -76,63 +76,43 @@ Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
 
 double log(double x)
 {
-	double hfsq,f,s,z,R,w,t1,t2,dk;
-	int32_t k,hx,i,j;
-	uint32_t lx;
-
-	EXTRACT_WORDS(hx, lx, x);
+	union {double f; uint64_t i;} u = {x};
+	double_t hfsq,f,s,z,R,w,t1,t2,dk;
+	uint32_t hx;
+	int k;
 
+	hx = u.i>>32;
 	k = 0;
-	if (hx < 0x00100000) {  /* x < 2**-1022  */
-		if (((hx&0x7fffffff)|lx) == 0)
-			return -two54/0.0;  /* log(+-0)=-inf */
-		if (hx < 0)
-			return (x-x)/0.0;   /* log(-#) = NaN */
-		/* subnormal number, scale up x */
+	if (hx < 0x00100000 || hx>>31) {
+		if (u.i<<1 == 0)
+			return -1/(x*x);  /* log(+-0)=-inf */
+		if (hx>>31)
+			return (x-x)/0.0; /* log(-#) = NaN */
+		/* subnormal number, scale x up */
 		k -= 54;
-		x *= two54;
-		GET_HIGH_WORD(hx,x);
-	}
-	if (hx >= 0x7ff00000)
-		return x+x;
-	k += (hx>>20) - 1023;
-	hx &= 0x000fffff;
-	i = (hx+0x95f64)&0x100000;
-	SET_HIGH_WORD(x, hx|(i^0x3ff00000));  /* normalize x or x/2 */
-	k += i>>20;
+		x *= 0x1p54;
+		u.f = x;
+		hx = u.i>>32;
+	} else if (hx >= 0x7ff00000) {
+		return x;
+	} else if (hx == 0x3ff00000 && u.i<<32 == 0)
+		return 0;
+
+	/* reduce x into [sqrt(2)/2, sqrt(2)] */
+	hx += 0x3ff00000 - 0x3fe6a09e;
+	k += (int)(hx>>20) - 0x3ff;
+	hx = (hx&0x000fffff) + 0x3fe6a09e;
+	u.i = (uint64_t)hx<<32 | (u.i&0xffffffff);
+	x = u.f;
+
 	f = x - 1.0;
-	if ((0x000fffff&(2+hx)) < 3) {  /* -2**-20 <= f < 2**-20 */
-		if (f == 0.0) {
-			if (k == 0) {
-				return 0.0;
-			}
-			dk = (double)k;
-			return dk*ln2_hi + dk*ln2_lo;
-		}
-		R = f*f*(0.5-0.33333333333333333*f);
-		if (k == 0)
-			return f - R;
-		dk = (double)k;
-		return dk*ln2_hi - ((R-dk*ln2_lo)-f);
-	}
+	hfsq = 0.5*f*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));
-		return dk*ln2_hi - ((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
-	} else {
-		if (k == 0)
-			return f - s*(f-R);
-		return dk*ln2_hi - ((s*(f-R)-dk*ln2_lo)-f);
-	}
+	dk = k;
+	return s*(hfsq+R) + dk*ln2_lo - hfsq + f + dk*ln2_hi;
 }