about summary refs log tree commit diff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/math/expl.c8
-rw-r--r--src/math/expm1l.c6
-rw-r--r--src/math/fma.c12
-rw-r--r--src/math/fmal.c12
-rw-r--r--src/math/log10l.c22
-rw-r--r--src/math/log2l.c20
-rw-r--r--src/math/logl.c20
-rw-r--r--src/math/powl.c103
8 files changed, 102 insertions, 101 deletions
diff --git a/src/math/expl.c b/src/math/expl.c
index 9507fd2e..b289ffec 100644
--- a/src/math/expl.c
+++ b/src/math/expl.c
@@ -102,13 +102,13 @@ long double expl(long double x)
 	if (x > MAXLOGL)
 		return INFINITY;
 	if (x < MINLOGL)
-		return 0.0L;
+		return 0.0;
 
 	/* Express e**x = e**g 2**n
 	 *   = e**g e**(n loge(2))
 	 *   = e**(g + n loge(2))
 	 */
-	px = floorl(LOG2EL * x + 0.5L); /* floor() truncates toward -infinity. */
+	px = floorl(LOG2EL * x + 0.5); /* floor() truncates toward -infinity. */
 	n = px;
 	x -= px * C1;
 	x -= px * C2;
@@ -120,8 +120,8 @@ long double expl(long double x)
 	xx = x * x;
 	px = x * __polevll(xx, P, 2);
 	x =  px/(__polevll(xx, Q, 3) - px);
-	x = 1.0L + ldexpl(x, 1);
-	x = ldexpl(x, n);
+	x = 1.0 + 2.0 * x;
+	x = scalbnl(x, n);
 	return x;
 }
 #endif
diff --git a/src/math/expm1l.c b/src/math/expm1l.c
index 2f94dfa2..ca0d141f 100644
--- a/src/math/expm1l.c
+++ b/src/math/expm1l.c
@@ -97,11 +97,11 @@ long double expm1l(long double x)
 		return x;
 	/* Minimum value.*/
 	if (x < minarg)
-		return -1.0L;
+		return -1.0;
 
 	xx = C1 + C2;
 	/* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
-	px = floorl (0.5 + x / xx);
+	px = floorl(0.5 + x / xx);
 	k = px;
 	/* remainder times ln 2 */
 	x -= px * C1;
@@ -116,7 +116,7 @@ long double expm1l(long double x)
 	/* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
 	 We have qx = exp(remainder ln 2) - 1, so
 	 exp(x) - 1  =  2^k (qx + 1) - 1  =  2^k qx + 2^k - 1.  */
-	px = ldexpl(1.0L, k);
+	px = scalbnl(1.0, k);
 	x = px * qx + (px - 1.0);
 	return x;
 }
diff --git a/src/math/fma.c b/src/math/fma.c
index 87d450c7..5fb95406 100644
--- a/src/math/fma.c
+++ b/src/math/fma.c
@@ -247,7 +247,7 @@ static inline double add_and_denormalize(double a, double b, int scale)
 			INSERT_WORD64(sum.hi, hibits);
 		}
 	}
-	return (ldexp(sum.hi, scale));
+	return scalbn(sum.hi, scale);
 }
 
 /*
@@ -364,7 +364,7 @@ double fma(double x, double y, double z)
 		}
 	}
 	if (spread <= DBL_MANT_DIG * 2)
-		zs = ldexp(zs, -spread);
+		zs = scalbn(zs, -spread);
 	else
 		zs = copysign(DBL_MIN, zs);
 
@@ -390,7 +390,7 @@ double fma(double x, double y, double z)
 		 */
 		fesetround(oround);
 		volatile double vzs = zs; /* XXX gcc CSE bug workaround */
-		return (xy.hi + vzs + ldexp(xy.lo, spread));
+		return xy.hi + vzs + scalbn(xy.lo, spread);
 	}
 
 	if (oround != FE_TONEAREST) {
@@ -400,13 +400,13 @@ double fma(double x, double y, double z)
 		 */
 		fesetround(oround);
 		adj = r.lo + xy.lo;
-		return (ldexp(r.hi + adj, spread));
+		return scalbn(r.hi + adj, spread);
 	}
 
 	adj = add_adjusted(r.lo, xy.lo);
 	if (spread + ilogb(r.hi) > -1023)
-		return (ldexp(r.hi + adj, spread));
+		return scalbn(r.hi + adj, spread);
 	else
-		return (add_and_denormalize(r.hi, adj, spread));
+		return add_and_denormalize(r.hi, adj, spread);
 }
 #endif
diff --git a/src/math/fmal.c b/src/math/fmal.c
index cbaf46eb..be64f145 100644
--- a/src/math/fmal.c
+++ b/src/math/fmal.c
@@ -115,7 +115,7 @@ static inline long double add_and_denormalize(long double a, long double b, int
 		if (bits_lost != 1 ^ (int)(u.bits.manl & 1))
 			sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
 	}
-	return (ldexp(sum.hi, scale));
+	return scalbnl(sum.hi, scale);
 }
 
 /*
@@ -228,7 +228,7 @@ long double fmal(long double x, long double y, long double z)
 		}
 	}
 	if (spread <= LDBL_MANT_DIG * 2)
-		zs = ldexpl(zs, -spread);
+		zs = scalbnl(zs, -spread);
 	else
 		zs = copysignl(LDBL_MIN, zs);
 
@@ -254,7 +254,7 @@ long double fmal(long double x, long double y, long double z)
 		 */
 		fesetround(oround);
 		volatile long double vzs = zs; /* XXX gcc CSE bug workaround */
-		return (xy.hi + vzs + ldexpl(xy.lo, spread));
+		return xy.hi + vzs + scalbnl(xy.lo, spread);
 	}
 
 	if (oround != FE_TONEAREST) {
@@ -264,13 +264,13 @@ long double fmal(long double x, long double y, long double z)
 		 */
 		fesetround(oround);
 		adj = r.lo + xy.lo;
-		return (ldexpl(r.hi + adj, spread));
+		return scalbnl(r.hi + adj, spread);
 	}
 
 	adj = add_adjusted(r.lo, xy.lo);
 	if (spread + ilogbl(r.hi) > -16383)
-		return (ldexpl(r.hi + adj, spread));
+		return scalbnl(r.hi + adj, spread);
 	else
-		return (add_and_denormalize(r.hi, adj, spread));
+		return add_and_denormalize(r.hi, adj, spread);
 }
 #endif
diff --git a/src/math/log10l.c b/src/math/log10l.c
index b954cc77..f0eeeafb 100644
--- a/src/math/log10l.c
+++ b/src/math/log10l.c
@@ -123,9 +123,9 @@ long double log10l(long double x)
 
 	if (isnan(x))
 		return x;
-	if(x <= 0.0L) {
-		if(x == 0.0L)
-			return -1.0L / (x - x);
+	if(x <= 0.0) {
+		if(x == 0.0)
+			return -1.0 / (x - x);
 		return (x - x) / (x - x);
 	}
 	if (x == INFINITY)
@@ -142,12 +142,12 @@ long double log10l(long double x)
 	if (e > 2 || e < -2) {
 		if (x < SQRTH) {  /* 2(2x-1)/(2x+1) */
 			e -= 1;
-			z = x - 0.5L;
-			y = 0.5L * z + 0.5L;
+			z = x - 0.5;
+			y = 0.5 * z + 0.5;
 		} else {  /*  2 (x-1)/(x+1)   */
-			z = x - 0.5L;
-			z -= 0.5L;
-			y = 0.5L * x  + 0.5L;
+			z = x - 0.5;
+			z -= 0.5;
+			y = 0.5 * x  + 0.5;
 		}
 		x = z / y;
 		z = x*x;
@@ -158,13 +158,13 @@ long double log10l(long double x)
 	/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 	if (x < SQRTH) {
 		e -= 1;
-		x = ldexpl(x, 1) - 1.0L; /*  2x - 1  */
+		x = 2.0*x - 1.0;
 	} else {
-		x = x - 1.0L;
+		x = x - 1.0;
 	}
 	z = x*x;
 	y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
-	y = y - ldexpl(z, -1);   /* -0.5x^2 + ... */
+	y = y - 0.5*z;
 
 done:
 	/* Multiply log of fraction by log10(e)
diff --git a/src/math/log2l.c b/src/math/log2l.c
index 4339c033..8ebce9c4 100644
--- a/src/math/log2l.c
+++ b/src/math/log2l.c
@@ -121,8 +121,8 @@ long double log2l(long double x)
 		return x;
 	if (x == INFINITY)
 		return x;
-	if (x <= 0.0L) {
-		if (x == 0.0L)
+	if (x <= 0.0) {
+		if (x == 0.0)
 			return -INFINITY;
 		return NAN;
 	}
@@ -139,12 +139,12 @@ long double log2l(long double x)
 	if (e > 2 || e < -2) {
 		if (x < SQRTH) {  /* 2(2x-1)/(2x+1) */
 			e -= 1;
-			z = x - 0.5L;
-			y = 0.5L * z + 0.5L;
+			z = x - 0.5;
+			y = 0.5 * z + 0.5;
 		} else {  /*  2 (x-1)/(x+1)   */
-			z = x - 0.5L;
-			z -= 0.5L;
-			y = 0.5L * x  + 0.5L;
+			z = x - 0.5;
+			z -= 0.5;
+			y = 0.5 * x + 0.5;
 		}
 		x = z / y;
 		z = x*x;
@@ -155,13 +155,13 @@ long double log2l(long double x)
 	/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 	if (x < SQRTH) {
 		e -= 1;
-		x = ldexpl(x, 1) - 1.0L; /*  2x - 1  */
+		x = 2.0*x - 1.0;
 	} else {
-		x = x - 1.0L;
+		x = x - 1.0;
 	}
 	z = x*x;
 	y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
-	y = y - ldexpl(z, -1);   /* -0.5x^2 + ... */
+	y = y - 0.5*z;
 
 done:
 	/* Multiply log of fraction by log2(e)
diff --git a/src/math/logl.c b/src/math/logl.c
index ee7ca64a..ffd81038 100644
--- a/src/math/logl.c
+++ b/src/math/logl.c
@@ -119,8 +119,8 @@ long double logl(long double x)
 		return x;
 	if (x == INFINITY)
 		return x;
-	if (x <= 0.0L) {
-		if (x == 0.0L)
+	if (x <= 0.0) {
+		if (x == 0.0)
 			return -INFINITY;
 		return NAN;
 	}
@@ -137,12 +137,12 @@ long double logl(long double x)
 	if (e > 2 || e < -2) {
 		if (x < SQRTH) {  /* 2(2x-1)/(2x+1) */
 			e -= 1;
-			z = x - 0.5L;
-			y = 0.5L * z + 0.5L;
+			z = x - 0.5;
+			y = 0.5 * z + 0.5;
 		} else {  /*  2 (x-1)/(x+1)   */
-			z = x - 0.5L;
-			z -= 0.5L;
-			y = 0.5L * x  + 0.5L;
+			z = x - 0.5;
+			z -= 0.5;
+			y = 0.5 * x  + 0.5;
 		}
 		x = z / y;
 		z = x*x;
@@ -156,14 +156,14 @@ long double logl(long double x)
 	/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
 	if (x < SQRTH) {
 		e -= 1;
-		x = ldexpl(x, 1) - 1.0L; /*  2x - 1  */
+		x = 2.0*x - 1.0;
 	} else {
-		x = x - 1.0L;
+		x = x - 1.0;
 	}
 	z = x*x;
 	y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6));
 	y = y + e * C2;
-	z = y - ldexpl(z, -1);   /*  y - 0.5 * z  */
+	z = y - 0.5*z;
 	/* Note, the sum of above terms does not exceed x/4,
 	 * so it contributes at most about 1/4 lsb to the error.
 	 */
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;
 }
generated by cgit-pink 1.4.1 (git 2.36.1) at 2024-05-16 13:39:02 +0000