RepoMirrors/musl
1
0
Fork 0
mirror of git://git.musl-libc.org/musl synced 2025-02-16 02:46:53 +00:00
Code Issues Packages Projects Releases Wiki Activity

use scalbn or *2.0 instead of ldexp, fix fmal

Browse Source 
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.
This commit is contained in:
nsz 2012-03-19 22:57:58 +01:00
parent 01fdfd491b
commit 2786c7d216
8 changed files with 102 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/math/expl.c
View File

@ -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

6
src/math/expm1l.c
View File

@ -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;
}

12
src/math/fma.c
View File

@ -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

12
src/math/fmal.c
View File

@ -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

22
src/math/log10l.c
View File

@ -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)

20
src/math/log2l.c
View File

@ -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)

20
src/math/logl.c
View File

@ -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.
*/

103
src/math/powl.c
View File

@ -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;
}