projects
/
musl
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
use alternate formula for acos asm to avoid loss of precision
[musl]
/
src
/
math
/
powl.c
diff --git
a/src/math/powl.c
b/src/math/powl.c
index
a6ee275
..
a1d2f07
100644
(file)
--- 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;
volatile long double z=0;
long double w=0, W=0, Wa=0, Wb=0, ya=0, yb=0, u=0;
- if (y == 0.0
L
)
- return 1.0
L
;
+ 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;
- if (y == 1.0
L
)
+ if (y == 1.0)
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.0
L || 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.0
L
)
- return 1.0
L
;
+ if (x == 1.0)
+ return 1.0;
if (y >= LDBL_MAX) {
if (y >= LDBL_MAX) {
- if (x > 1.0
L
)
+ if (x > 1.0)
return INFINITY;
return INFINITY;
- if (x > 0.0
L && x < 1.0L
)
- return 0.0
L
;
- if (x < -1.0
L
)
+ if (x > 0.0
&& x < 1.0
)
+ return 0.0;
+ if (x < -1.0)
return INFINITY;
return INFINITY;
- if (x > -1.0
L && x < 0.0L
)
- return 0.0
L
;
+ if (x > -1.0
&& x < 0.0
)
+ return 0.0;
}
if (y <= -LDBL_MAX) {
}
if (y <= -LDBL_MAX) {
- if (x > 1.0
L
)
- return 0.0
L
;
- if (x > 0.0
L && x < 1.0L
)
+ if (x > 1.0)
+ return 0.0;
+ if (x > 0.0
&& x < 1.0
)
return INFINITY;
return INFINITY;
- if (x < -1.0
L
)
- return 0.0
L
;
- if (x > -1.0
L && 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) {
- if (y > 0.0
L
)
+ if (y > 0.0)
return INFINITY;
return INFINITY;
- return 0.0
L
;
+ return 0.0;
}
w = floorl(y);
}
w = floorl(y);
@@
-253,29
+253,29
@@
long double powl(long double x, long double y)
yoddint = 0;
if (iyflg) {
ya = fabsl(y);
yoddint = 0;
if (iyflg) {
ya = fabsl(y);
- ya = floorl(0.5
L
* ya);
- yb = 0.5
L
* 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 (y > 0.0
L
) {
+ if (y > 0.0) {
if (yoddint)
return -INFINITY;
return INFINITY;
}
if (yoddint)
return -INFINITY;
return INFINITY;
}
- if (y < 0.0
L
) {
+ if (y < 0.0) {
if (yoddint)
if (yoddint)
- return -0.0
L
;
+ return -0.0;
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.0
L
) {
- if (x == 0.0
L
) {
+ 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;
@@
-283,12
+283,12
@@
long double powl(long double x, long double y)
}
if (y > 0.0) {
if (signbit(x) && yoddint)
}
if (y > 0.0) {
if (signbit(x) && yoddint)
- return -0.0
L
;
+ return -0.0;
return 0.0;
}
return 0.0;
}
- if (y == 0.0
L
)
- return 1.0
L
; /* 0**0 */
- return 0.0
L
; /* 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 */
@@
-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));
*/
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
/* 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;
/* 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. */
@@
-380,7
+381,7
@@
long double powl(long double x, long double y)
H = Fb + Gb;
Ha = reducl(H);
H = Fb + Gb;
Ha = reducl(H);
- w =
ldexp
l( Ga+Ha, LNXT );
+ w =
scalbn
l( Ga+Ha, LNXT );
/* Test the power of 2 for overflow */
if (w > MEXP)
/* 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;
e = w;
Hb = H - Ha;
- if (Hb > 0.0
L
) {
+ if (Hb > 0.0) {
e += 1;
e += 1;
- Hb -= 1.0
L
/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.
@@
-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;
w = douba(e);
z = w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */
z = z + w;
- z =
ldexp
l(z, i); /* multiply by integer power of 2 */
+ z =
scalbn
l(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.
*/
- w =
ldexpl(y, -1)
;
+ w =
0.5*y
;
w = floorl(w);
w = floorl(w);
- w =
ldexpl(w, 1)
;
+ w =
2.0*w
;
if (w != y)
z = -z; /* odd exponent */
}
if (w != y)
z = -z; /* odd exponent */
}
@@
-438,9
+439,9
@@
static long double reducl(long double x)
{
long double t;
{
long double t;
- t =
ldexp
l(x, LNXT);
+ t =
scalbn
l(x, LNXT);
t = floorl(t);
t = floorl(t);
- t =
ldexp
l(t, -LNXT);
+ t =
scalbn
l(t, -LNXT);
return t;
}
return t;
}
@@
-483,18
+484,18
@@
static long double powil(long double x, int nn)
long double s;
int n, e, sign, asign, lx;
long double s;
int n, e, sign, asign, lx;
- if (x == 0.0
L
) {
+ if (x == 0.0) {
if (nn == 0)
if (nn == 0)
- return 1.0
L
;
+ return 1.0;
else if (nn < 0)
return LDBL_MAX;
else if (nn < 0)
return LDBL_MAX;
- return 0.0
L
;
+ return 0.0;
}
if (nn == 0)
}
if (nn == 0)
- return 1.0
L
;
+ return 1.0;
- if (x < 0.0
L
) {
+ if (x < 0.0) {
asign = -1;
x = -x;
} else
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);
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.5
L
+ lx) * nn * LOGE2L;
+ s = (2.9142135623730950L * s - 0.5 + lx) * nn * LOGE2L;
} else {
s = LOGE2L * e;
}
} 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.
*/
* since roundoff error in 1.0/x will be amplified.
* The precise demarcation should be the gradual underflow threshold.
*/
- if (s < -MAXLOGL+2.0
L
) {
- x = 1.0
L
/x;
+ if (s < -MAXLOGL+2.0) {
+ x = 1.0/x;
sign = -sign;
}
sign = -sign;
}
@@
-539,7
+540,7
@@
static long double powil(long double x, int nn)
if (n & 1)
y = x;
else {
if (n & 1)
y = x;
else {
- y = 1.0
L
;
+ y = 1.0;
asign = 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)
if (asign)
y = -y; /* odd power of negative number */
if (sign < 0)
- y = 1.0
L
/y;
+ y = 1.0/y;
return y;
}
return y;
}