X-Git-Url: http://nsz.repo.hu/git/?p=musl;a=blobdiff_plain;f=src%2Fmath%2Fpowl.c;h=ce6274cf73e00c5233acd4954e3d7d2dccf6ff70;hp=690f2942a034ed02dad81dd92567bdde1d9ed446;hb=04ccbdca6d88738e23e0d6a622ad33854c468646;hpb=b69f695acedd4ce2798ef9ea28d834ceccc789bd

diff --git a/src/math/powl.c b/src/math/powl.c
index 690f2942..ce6274cf 100644
--- a/src/math/powl.c
+++ b/src/math/powl.c
@@ -78,19 +78,17 @@ long double powl(long double x, long double y)
 
 /* Table size */
 #define NXT 32
-/* log2(Table size) */
-#define LNXT 5
 
 /* log(1+x) =  x - .5x^2 + x^3 *  P(z)/Q(z)
  * on the domain  2^(-1/32) - 1  <=  x  <=  2^(1/32) - 1
  */
-static long double P[] = {
+static const long double P[] = {
  8.3319510773868690346226E-4L,
  4.9000050881978028599627E-1L,
  1.7500123722550302671919E0L,
  1.4000100839971580279335E0L,
 };
-static long double Q[] = {
+static const long double Q[] = {
 /* 1.0000000000000000000000E0L,*/
  5.2500282295834889175431E0L,
  8.4000598057587009834666E0L,
@@ -99,7 +97,7 @@ static long double Q[] = {
 /* A[i] = 2^(-i/32), rounded to IEEE long double precision.
  * If i is even, A[i] + B[i/2] gives additional accuracy.
  */
-static long double A[33] = {
+static const long double A[33] = {
  1.0000000000000000000000E0L,
  9.7857206208770013448287E-1L,
  9.5760328069857364691013E-1L,
@@ -134,7 +132,7 @@ static long double A[33] = {
  5.1094857432705833910408E-1L,
  5.0000000000000000000000E-1L,
 };
-static long double B[17] = {
+static const long double B[17] = {
  0.0000000000000000000000E0L,
  2.6176170809902549338711E-20L,
 -1.0126791927256478897086E-20L,
@@ -157,7 +155,7 @@ static long double B[17] = {
 /* 2^x = 1 + x P(x),
  * on the interval -1/32 <= x <= 0
  */
-static long double R[] = {
+static const long double R[] = {
  1.5089970579127659901157E-5L,
  1.5402715328927013076125E-4L,
  1.3333556028915671091390E-3L,
@@ -167,8 +165,6 @@ static long double R[] = {
  6.9314718055994530931447E-1L,
 };
 
-#define douba(k) A[k]
-#define doubb(k) B[k]
 #define MEXP (NXT*16384.0L)
 /* The following if denormal numbers are supported, else -MEXP: */
 #define MNEXP (-NXT*(16384.0L+64.0L))
@@ -188,11 +184,9 @@ static long double R[] = {
 static const long double MAXLOGL = 1.1356523406294143949492E4L;
 static const long double MINLOGL = -1.13994985314888605586758E4L;
 static const long double LOGE2L = 6.9314718055994530941723E-1L;
-static volatile long double z;
-static long double w, W, Wa, Wb, ya, yb, u;
 static const long double huge = 0x1p10000L;
 /* XXX Prevent gcc from erroneously constant folding this. */
-static volatile long double twom10000 = 0x1p-10000L;
+static const volatile long double twom10000 = 0x1p-10000L;
 
 static long double reducl(long double);
 static long double powil(long double, int);
@@ -202,48 +196,48 @@ long double powl(long double x, long double y)
 	/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
 	int i, nflg, iyflg, yoddint;
 	long e;
+	volatile long double z=0;
+	long double w=0, W=0, Wa=0, Wb=0, ya=0, yb=0, u=0;
 
-	if (y == 0.0L)
-		return 1.0L;
-	if (isnan(x))
+	/* make sure no invalid exception is raised by nan comparision */
+	if (isnan(x)) {
+		if (!isnan(y) && y == 0.0)
+			return 1.0;
 		return x;
-	if (isnan(y))
+	}
+	if (isnan(y)) {
+		if (x == 1.0)
+			return 1.0;
 		return y;
-	if (y == 1.0L)
+	}
+	if (x == 1.0)
+		return 1.0; /* 1**y = 1, even if y is nan */
+	if (x == -1.0 && !isfinite(y))
+		return 1.0; /* -1**inf = 1 */
+	if (y == 0.0)
+		return 1.0; /* x**0 = 1, even if x is nan */
+	if (y == 1.0)
 		return x;
-
-	// FIXME: this is wrong, see pow special cases in c99 F.9.4.4
-	if (!isfinite(y) && (x == -1.0L || x == 1.0L) )
-		return y - y;   /* +-1**inf is NaN */
-	if (x == 1.0L)
-		return 1.0L;
 	if (y >= LDBL_MAX) {
-		if (x > 1.0L)
-			return INFINITY;
-		if (x > 0.0L && x < 1.0L)
-			return 0.0L;
-		if (x < -1.0L)
+		if (x > 1.0 || x < -1.0)
 			return INFINITY;
-		if (x > -1.0L && x < 0.0L)
-			return 0.0L;
+		if (x != 0.0)
+			return 0.0;
 	}
 	if (y <= -LDBL_MAX) {
-		if (x > 1.0L)
-			return 0.0L;
-		if (x > 0.0L && x < 1.0L)
-			return INFINITY;
-		if (x < -1.0L)
-			return 0.0L;
-		if (x > -1.0L && x < 0.0L)
+		if (x > 1.0 || x < -1.0)
+			return 0.0;
+		if (x != 0.0)
 			return INFINITY;
 	}
 	if (x >= LDBL_MAX) {
-		if (y > 0.0L)
+		if (y > 0.0)
 			return INFINITY;
-		return 0.0L;
+		return 0.0;
 	}
 
 	w = floorl(y);
+
 	/* Set iyflg to 1 if y is an integer. */
 	iyflg = 0;
 	if (w == y)
@@ -253,76 +247,66 @@ long double powl(long double x, long double y)
 	yoddint = 0;
 	if (iyflg) {
 		ya = fabsl(y);
-		ya = floorl(0.5L * ya);
-		yb = 0.5L * fabsl(w);
+		ya = floorl(0.5 * ya);
+		yb = 0.5 * fabsl(w);
 		if( ya != yb )
 			yoddint = 1;
 	}
 
 	if (x <= -LDBL_MAX) {
-		if (y > 0.0L) {
+		if (y > 0.0) {
 			if (yoddint)
 				return -INFINITY;
 			return INFINITY;
 		}
-		if (y < 0.0L) {
+		if (y < 0.0) {
 			if (yoddint)
-				return -0.0L;
+				return -0.0;
 			return 0.0;
 		}
 	}
-
-
-	nflg = 0;       /* flag = 1 if x<0 raised to integer power */
-	if (x <= 0.0L) {
-		if (x == 0.0L) {
+	nflg = 0; /* (x<0)**(odd int) */
+	if (x <= 0.0) {
+		if (x == 0.0) {
 			if (y < 0.0) {
 				if (signbit(x) && yoddint)
-					return -INFINITY;
-				return INFINITY;
+					/* (-0.0)**(-odd int) = -inf, divbyzero */
+					return -1.0/0.0;
+				/* (+-0.0)**(negative) = inf, divbyzero */
+				return 1.0/0.0;
 			}
-			if (y > 0.0) {
-				if (signbit(x) && yoddint)
-					return -0.0L;
-				return 0.0;
-			}
-			if (y == 0.0L)
-				return 1.0L;  /*   0**0   */
-			return 0.0L;  /*   0**y   */
+			if (signbit(x) && yoddint)
+				return -0.0;
+			return 0.0;
 		}
 		if (iyflg == 0)
 			return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */
-		nflg = 1;
+		/* (x<0)**(integer) */
+		if (yoddint)
+			nflg = 1; /* negate result */
+		x = -x;
 	}
-
-	/* Integer power of an integer.  */
-	if (iyflg) {
-		i = w;
-		w = floorl(x);
-		if (w == x && fabsl(y) < 32768.0) {
-			w = powil(x, (int)y);
-			return w;
-		}
+	/* (+integer)**(integer)  */
+	if (iyflg && floorl(x) == x && fabsl(y) < 32768.0) {
+		w = powil(x, (int)y);
+		return nflg ? -w : w;
 	}
 
-	if (nflg)
-		x = fabsl(x);
-
 	/* separate significand from exponent */
 	x = frexpl(x, &i);
 	e = i;
 
 	/* find significand in antilog table A[] */
 	i = 1;
-	if (x <= douba(17))
+	if (x <= A[17])
 		i = 17;
-	if (x <= douba(i+8))
+	if (x <= A[i+8])
 		i += 8;
-	if (x <= douba(i+4))
+	if (x <= A[i+4])
 		i += 4;
-	if (x <= douba(i+2))
+	if (x <= A[i+2])
 		i += 2;
-	if (x >= douba(1))
+	if (x >= A[1])
 		i = -1;
 	i += 1;
 
@@ -333,9 +317,9 @@ long double powl(long double x, long double y)
 	 *
 	 * log(x/a) = log(1+v),  v = x/a - 1 = (x-a)/a
 	 */
-	x -= douba(i);
-	x -= doubb(i/2);
-	x /= douba(i);
+	x -= A[i];
+	x -= B[i/2];
+	x /= A[i];
 
 	/* rational approximation for log(1+v):
 	 *
@@ -343,7 +327,7 @@ long double powl(long double x, long double y)
 	 */
 	z = x*x;
 	w = x * (z * __polevll(x, P, 3) / __p1evll(x, Q, 3));
-	w = w - ldexpl(z, -1);  /*  w - 0.5 * z  */
+	w = w - 0.5*z;
 
 	/* Convert to base 2 logarithm:
 	 * multiply by log2(e) = 1 + LOG2EA
@@ -355,7 +339,7 @@ long double powl(long double x, long double y)
 
 	/* Compute exponent term of the base 2 logarithm. */
 	w = -i;
-	w = ldexpl(w, -LNXT); /* divide by NXT */
+	w /= NXT;
 	w += e;
 	/* Now base 2 log of x is w + z. */
 
@@ -380,7 +364,7 @@ long double powl(long double x, long double y)
 
 	H = Fb + Gb;
 	Ha = reducl(H);
-	w = ldexpl( Ga+Ha, LNXT );
+	w = (Ga + Ha) * NXT;
 
 	/* Test the power of 2 for overflow */
 	if (w > MEXP)
@@ -391,9 +375,9 @@ long double powl(long double x, long double y)
 	e = w;
 	Hb = H - Ha;
 
-	if (Hb > 0.0L) {
+	if (Hb > 0.0) {
 		e += 1;
-		Hb -= 1.0L/NXT;  /*0.0625L;*/
+		Hb -= 1.0/NXT;  /*0.0625L;*/
 	}
 
 	/* Now the product y * log2(x)  =  Hb + e/NXT.
@@ -412,23 +396,13 @@ long double powl(long double x, long double y)
 		i = 1;
 	i = e/NXT + i;
 	e = NXT*i - e;
-	w = douba(e);
+	w = A[e];
 	z = w * z;  /*  2**-e * ( 1 + (2**Hb-1) )  */
 	z = z + w;
-	z = ldexpl(z, i);  /* multiply by integer power of 2 */
-
-	if (nflg) {
-		/* For negative x,
-		 * find out if the integer exponent
-		 * is odd or even.
-		 */
-		w = ldexpl(y, -1);
-		w = floorl(w);
-		w = ldexpl(w, 1);
-		if (w != y)
-			z = -z;  /* odd exponent */
-	}
+	z = scalbnl(z, i);  /* multiply by integer power of 2 */
 
+	if (nflg)
+		z = -z;
 	return z;
 }
 
@@ -438,15 +412,14 @@ static long double reducl(long double x)
 {
 	long double t;
 
-	t = ldexpl(x, LNXT);
+	t = x * NXT;
 	t = floorl(t);
-	t = ldexpl(t, -LNXT);
+	t = t / NXT;
 	return t;
 }
 
-/*                                                      powil.c
- *
- *      Real raised to integer power, long double precision
+/*
+ *      Positive real raised to integer power, long double precision
  *
  *
  * SYNOPSIS:
@@ -459,7 +432,7 @@ static long double reducl(long double x)
  *
  * DESCRIPTION:
  *
- * Returns argument x raised to the nth power.
+ * Returns argument x>0 raised to the nth power.
  * The routine efficiently decomposes n as a sum of powers of
  * two. The desired power is a product of two-to-the-kth
  * powers of x.  Thus to compute the 32767 power of x requires
@@ -481,24 +454,10 @@ static long double powil(long double x, int nn)
 {
 	long double ww, y;
 	long double s;
-	int n, e, sign, asign, lx;
-
-	if (x == 0.0L) {
-		if (nn == 0)
-			return 1.0L;
-		else if (nn < 0)
-			return LDBL_MAX;
-		return 0.0L;
-	}
+	int n, e, sign, lx;
 
 	if (nn == 0)
-		return 1.0L;
-
-	if (x < 0.0L) {
-		asign = -1;
-		x = -x;
-	} else
-		asign = 0;
+		return 1.0;
 
 	if (nn < 0) {
 		sign = -1;
@@ -516,7 +475,7 @@ static long double powil(long double x, int nn)
 	e = (lx - 1)*n;
 	if ((e == 0) || (e > 64) || (e < -64)) {
 		s = (s - 7.0710678118654752e-1L) / (s +  7.0710678118654752e-1L);
-		s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L;
+		s = (2.9142135623730950L * s - 0.5 + lx) * nn * LOGE2L;
 	} else {
 		s = LOGE2L * e;
 	}
@@ -530,18 +489,16 @@ static long double powil(long double x, int nn)
 	 * since roundoff error in 1.0/x will be amplified.
 	 * The precise demarcation should be the gradual underflow threshold.
 	 */
-	if (s < -MAXLOGL+2.0L) {
-		x = 1.0L/x;
+	if (s < -MAXLOGL+2.0) {
+		x = 1.0/x;
 		sign = -sign;
 	}
 
 	/* First bit of the power */
 	if (n & 1)
 		y = x;
-	else {
-		y = 1.0L;
-		asign = 0;
-	}
+	else
+		y = 1.0;
 
 	ww = x;
 	n >>= 1;
@@ -552,10 +509,8 @@ static long double powil(long double x, int nn)
 		n >>= 1;
 	}
 
-	if (asign)
-		y = -y;  /* odd power of negative number */
 	if (sign < 0)
-		y = 1.0L/y;
+		y = 1.0/y;
 	return y;
 }