summaryrefslogtreecommitdiff
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 v1.2.1 (git 2.18.0) at 2024-04-26 15:28:35 +0000