summaryrefslogtreecommitdiff
path: root/src/math/powl.c
diff options
context:
space:
mode:
authornsz <nsz@port70.net>2012-03-19 22:57:58 +0100
committernsz <nsz@port70.net>2012-03-19 22:57:58 +0100
commit2786c7d21611b9fa3b2fe356542cf213e7dd0ba4 (patch)
treeb3954e9cec7580f5dc851491d3b60d808aae4259 /src/math/powl.c
parent01fdfd491b5d83b72099fbae14c4a71ed8e0b945 (diff)
downloadmusl-2786c7d21611b9fa3b2fe356542cf213e7dd0ba4.tar.gz
use scalbn or *2.0 instead of ldexp, fix fmal
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.
Diffstat (limited to 'src/math/powl.c')
-rw-r--r--src/math/powl.c103
1 files changed, 52 insertions, 51 deletions
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;
}