X-Git-Url: http://nsz.repo.hu/git/?a=blobdiff_plain;f=src%2Fmath%2Fpowl.c;h=a1d2f0769eb550ac39f1cb36f93903c55ebf8e10;hb=634c3a63027aa4a693b64fae0e2f6e1635558e93;hp=a6ee275f25eedb4c95f367b48ce893009702f8bc;hpb=9e2a895aaaa4a3985e94ae4f3e24c1af65f9bb34;p=musl

diff --git a/src/math/powl.c b/src/math/powl.c
index a6ee275f..a1d2f076 100644
--- a/src/math/powl.c
+++ b/src/math/powl.c
@@ -203,44 +203,44 @@ long double powl(long double x, long double y)
 	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 (y == 0.0)
+		return 1.0;
 	if (isnan(x))
 		return x;
 	if (isnan(y))
 		return y;
-	if (y == 1.0L)
+	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) )
+	if (!isfinite(y) && (x == -1.0 || x == 1.0) )
 		return y - y;   /* +-1**inf is NaN */
-	if (x == 1.0L)
-		return 1.0L;
+	if (x == 1.0)
+		return 1.0;
 	if (y >= LDBL_MAX) {
-		if (x > 1.0L)
+		if (x > 1.0)
 			return INFINITY;
-		if (x > 0.0L && x < 1.0L)
-			return 0.0L;
-		if (x < -1.0L)
+		if (x > 0.0 && x < 1.0)
+			return 0.0;
+		if (x < -1.0)
 			return INFINITY;
-		if (x > -1.0L && x < 0.0L)
-			return 0.0L;
+		if (x > -1.0 && x < 0.0)
+			return 0.0;
 	}
 	if (y <= -LDBL_MAX) {
-		if (x > 1.0L)
-			return 0.0L;
-		if (x > 0.0L && x < 1.0L)
+		if (x > 1.0)
+			return 0.0;
+		if (x > 0.0 && x < 1.0)
 			return INFINITY;
-		if (x < -1.0L)
-			return 0.0L;
-		if (x > -1.0L && x < 0.0L)
+		if (x < -1.0)
+			return 0.0;
+		if (x > -1.0 && 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);
@@ -253,29 +253,29 @@ 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) {
+	if (x <= 0.0) {
+		if (x == 0.0) {
 			if (y < 0.0) {
 				if (signbit(x) && yoddint)
 					return -INFINITY;
@@ -283,12 +283,12 @@ long double powl(long double x, long double y)
 			}
 			if (y > 0.0) {
 				if (signbit(x) && yoddint)
-					return -0.0L;
+					return -0.0;
 				return 0.0;
 			}
-			if (y == 0.0L)
-				return 1.0L;  /*   0**0   */
-			return 0.0L;  /*   0**y   */
+			if (y == 0.0)
+				return 1.0;  /*   0**0   */
+			return 0.0;  /*   0**y   */
 		}
 		if (iyflg == 0)
 			return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */
@@ -343,7 +343,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 +355,8 @@ 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 */
+	// TODO: use w * 0x1p-5;
+	w = scalbnl(w, -LNXT); /* divide by NXT */
 	w += e;
 	/* Now base 2 log of x is w + z. */
 
@@ -380,7 +381,7 @@ long double powl(long double x, long double y)
 
 	H = Fb + Gb;
 	Ha = reducl(H);
-	w = ldexpl( Ga+Ha, LNXT );
+	w = scalbnl( Ga+Ha, LNXT );
 
 	/* Test the power of 2 for overflow */
 	if (w > MEXP)
@@ -391,9 +392,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.
@@ -415,16 +416,16 @@ long double powl(long double x, long double y)
 	w = douba(e);
 	z = w * z;  /*  2**-e * ( 1 + (2**Hb-1) )  */
 	z = z + w;
-	z = ldexpl(z, i);  /* multiply by integer power of 2 */
+	z = scalbnl(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 = 0.5*y;
 		w = floorl(w);
-		w = ldexpl(w, 1);
+		w = 2.0*w;
 		if (w != y)
 			z = -z;  /* odd exponent */
 	}
@@ -438,9 +439,9 @@ static long double reducl(long double x)
 {
 	long double t;
 
-	t = ldexpl(x, LNXT);
+	t = scalbnl(x, LNXT);
 	t = floorl(t);
-	t = ldexpl(t, -LNXT);
+	t = scalbnl(t, -LNXT);
 	return t;
 }
 
@@ -483,18 +484,18 @@ static long double powil(long double x, int nn)
 	long double s;
 	int n, e, sign, asign, lx;
 
-	if (x == 0.0L) {
+	if (x == 0.0) {
 		if (nn == 0)
-			return 1.0L;
+			return 1.0;
 		else if (nn < 0)
 			return LDBL_MAX;
-		return 0.0L;
+		return 0.0;
 	}
 
 	if (nn == 0)
-		return 1.0L;
+		return 1.0;
 
-	if (x < 0.0L) {
+	if (x < 0.0) {
 		asign = -1;
 		x = -x;
 	} else
@@ -516,7 +517,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,8 +531,8 @@ 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;
 	}
 
@@ -539,7 +540,7 @@ static long double powil(long double x, int nn)
 	if (n & 1)
 		y = x;
 	else {
-		y = 1.0L;
+		y = 1.0;
 		asign = 0;
 	}
 
@@ -555,7 +556,7 @@ static long double powil(long double x, int nn)
 	if (asign)
 		y = -y;  /* odd power of negative number */
 	if (sign < 0)
-		y = 1.0L/y;
+		y = 1.0/y;
 	return y;
 }