Some code assumed ldexp(x, 1) is faster than 2.0*x,
but ldexp is a wrapper around scalbn which uses
multiplications inside, so this optimization is
wrong.
This commit also fixes fmal which accidentally
used ldexp instead of ldexpl loosing precision.
There are various additional changes from the
work-in-progress const cleanups.
if (x > MAXLOGL)
return INFINITY;
if (x < MINLOGL)
if (x > MAXLOGL)
return INFINITY;
if (x < MINLOGL)
/* Express e**x = e**g 2**n
* = e**g e**(n loge(2))
* = e**(g + n loge(2))
*/
/* 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;
n = px;
x -= px * C1;
x -= px * C2;
xx = x * x;
px = x * __polevll(xx, P, 2);
x = px/(__polevll(xx, Q, 3) - px);
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;
/* Minimum value.*/
if (x < minarg)
return x;
/* Minimum value.*/
if (x < minarg)
xx = C1 + C2;
/* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
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;
k = px;
/* remainder times ln 2 */
x -= px * C1;
/* 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. */
/* 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. */
x = px * qx + (px - 1.0);
return x;
}
x = px * qx + (px - 1.0);
return x;
}
INSERT_WORD64(sum.hi, hibits);
}
}
INSERT_WORD64(sum.hi, hibits);
}
}
- return (ldexp(sum.hi, scale));
+ return scalbn(sum.hi, scale);
}
}
if (spread <= DBL_MANT_DIG * 2)
}
}
if (spread <= DBL_MANT_DIG * 2)
- zs = ldexp(zs, -spread);
+ zs = scalbn(zs, -spread);
else
zs = copysign(DBL_MIN, zs);
else
zs = copysign(DBL_MIN, zs);
*/
fesetround(oround);
volatile double vzs = zs; /* XXX gcc CSE bug workaround */
*/
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) {
}
if (oround != FE_TONEAREST) {
*/
fesetround(oround);
adj = r.lo + xy.lo;
*/
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)
}
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);
- return (add_and_denormalize(r.hi, adj, spread));
+ return add_and_denormalize(r.hi, adj, spread);
if (bits_lost != 1 ^ (int)(u.bits.manl & 1))
sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
}
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);
}
}
if (spread <= LDBL_MANT_DIG * 2)
}
}
if (spread <= LDBL_MANT_DIG * 2)
- zs = ldexpl(zs, -spread);
+ zs = scalbnl(zs, -spread);
else
zs = copysignl(LDBL_MIN, zs);
else
zs = copysignl(LDBL_MIN, zs);
*/
fesetround(oround);
volatile long double vzs = zs; /* XXX gcc CSE bug workaround */
*/
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) {
}
if (oround != FE_TONEAREST) {
*/
fesetround(oround);
adj = r.lo + xy.lo;
*/
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)
}
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);
- return (add_and_denormalize(r.hi, adj, spread));
+ return add_and_denormalize(r.hi, adj, spread);
- 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)
return (x - x) / (x - x);
}
if (x == INFINITY)
if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1;
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) */
} 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;
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if (x < SQRTH) {
e -= 1;
/* 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 */
}
z = x*x;
y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
}
z = x*x;
y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
- y = y - ldexpl(z, -1); /* -0.5x^2 + ... */
done:
/* Multiply log of fraction by log10(e)
done:
/* Multiply log of fraction by log10(e)
return x;
if (x == INFINITY)
return 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;
}
return -INFINITY;
return NAN;
}
if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1;
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) */
} 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;
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if (x < SQRTH) {
e -= 1;
/* 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 */
}
z = x*x;
y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
}
z = x*x;
y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 7));
- y = y - ldexpl(z, -1); /* -0.5x^2 + ... */
done:
/* Multiply log of fraction by log2(e)
done:
/* Multiply log of fraction by log2(e)
return x;
if (x == INFINITY)
return 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;
}
return -INFINITY;
return NAN;
}
if (e > 2 || e < -2) {
if (x < SQRTH) { /* 2(2x-1)/(2x+1) */
e -= 1;
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) */
} 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;
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if (x < SQRTH) {
e -= 1;
/* 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 */
}
z = x*x;
y = x * (z * __polevll(x, P, 6) / __p1evll(x, Q, 6));
y = y + e * C2;
}
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 */
/* Note, the sum of above terms does not exceed x/4,
* so it contributes at most about 1/4 lsb to the error.
*/
/* Note, the sum of above terms does not exceed x/4,
* so it contributes at most about 1/4 lsb to the error.
*/
volatile long double z=0;
long double w=0, W=0, Wa=0, Wb=0, ya=0, yb=0, u=0;
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 (isnan(x))
return x;
if (isnan(y))
return y;
return x;
// FIXME: this is wrong, see pow special cases in c99 F.9.4.4
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 */
return y - y; /* +-1**inf is NaN */
- if (x == 1.0L)
- return 1.0L;
+ if (x == 1.0)
+ return 1.0;
- 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)
- if (x > -1.0L && x < 0.0L)
- return 0.0L;
+ if (x > -1.0 && x < 0.0)
+ return 0.0;
- 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)
- 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) {
return INFINITY;
}
if (x >= LDBL_MAX) {
yoddint = 0;
if (iyflg) {
ya = fabsl(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( ya != yb )
yoddint = 1;
}
if (x <= -LDBL_MAX) {
if (yoddint)
return -INFINITY;
return INFINITY;
}
if (yoddint)
return -INFINITY;
return INFINITY;
}
return 0.0;
}
}
nflg = 0; /* flag = 1 if x<0 raised to integer power */
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;
if (y < 0.0) {
if (signbit(x) && yoddint)
return -INFINITY;
}
if (y > 0.0) {
if (signbit(x) && yoddint)
}
if (y > 0.0) {
if (signbit(x) && yoddint)
- 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 */
}
if (iyflg == 0)
return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */
*/
z = x*x;
w = x * (z * __polevll(x, P, 3) / __p1evll(x, Q, 3));
*/
z = x*x;
w = x * (z * __polevll(x, P, 3) / __p1evll(x, Q, 3));
- w = w - ldexpl(z, -1); /* w - 0.5 * z */
/* Convert to base 2 logarithm:
* multiply by log2(e) = 1 + LOG2EA
/* Convert to base 2 logarithm:
* multiply by log2(e) = 1 + LOG2EA
/* Compute exponent term of the base 2 logarithm. */
w = -i;
/* 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. */
w += e;
/* Now base 2 log of x is w + z. */
H = Fb + Gb;
Ha = reducl(H);
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)
/* Test the power of 2 for overflow */
if (w > MEXP)
- Hb -= 1.0L/NXT; /*0.0625L;*/
+ Hb -= 1.0/NXT; /*0.0625L;*/
}
/* Now the product y * log2(x) = Hb + e/NXT.
}
/* Now the product y * log2(x) = Hb + e/NXT.
w = douba(e);
z = w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */
z = z + w;
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.
*/
if (nflg) {
/* For negative x,
* find out if the integer exponent
* is odd or even.
*/
if (w != y)
z = -z; /* odd exponent */
}
if (w != y)
z = -z; /* odd exponent */
}
long double s;
int n, e, sign, asign, lx;
long double s;
int n, e, sign, asign, lx;
else if (nn < 0)
return LDBL_MAX;
else if (nn < 0)
return LDBL_MAX;
asign = -1;
x = -x;
} else
asign = -1;
x = -x;
} else
e = (lx - 1)*n;
if ((e == 0) || (e > 64) || (e < -64)) {
s = (s - 7.0710678118654752e-1L) / (s + 7.0710678118654752e-1L);
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;
}
} else {
s = LOGE2L * e;
}
* since roundoff error in 1.0/x will be amplified.
* The precise demarcation should be the gradual underflow threshold.
*/
* 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;
if (asign)
y = -y; /* odd power of negative number */
if (sign < 0)
if (asign)
y = -y; /* odd power of negative number */
if (sign < 0)