projects
/
musl
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
move register_t and u_int64_t (back) to alltypes
[musl]
/
src
/
math
/
fmal.c
diff --git
a/src/math/fmal.c
b/src/math/fmal.c
index
200bd5a
..
87e30fc
100644
(file)
--- 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);
}
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
);
}
/*
}
/*
@@
-162,6
+162,7
@@
static inline struct dd dd_mul(long double a, long double b)
*/
long double fmal(long double x, long double y, long double z)
{
*/
long double fmal(long double x, long double y, long double z)
{
+ #pragma STDC FENV_ACCESS ON
long double xs, ys, zs, adj;
struct dd xy, r;
int oround;
long double xs, ys, zs, adj;
struct dd xy, r;
int oround;
@@
-173,14
+174,14
@@
long double fmal(long double x, long double y, long double z)
* return values here are crucial in handling special cases involving
* infinities, NaNs, overflows, and signed zeroes correctly.
*/
* return values here are crucial in handling special cases involving
* infinities, NaNs, overflows, and signed zeroes correctly.
*/
- if (x == 0.0 || y == 0.0)
- return (x * y + z);
- if (z == 0.0)
- return (x * y);
if (!isfinite(x) || !isfinite(y))
return (x * y + z);
if (!isfinite(z))
return (z);
if (!isfinite(x) || !isfinite(y))
return (x * y + z);
if (!isfinite(z))
return (z);
+ if (x == 0.0 || y == 0.0)
+ return (x * y + z);
+ if (z == 0.0)
+ return (x * y);
xs = frexpl(x, &ex);
ys = frexpl(y, &ey);
xs = frexpl(x, &ex);
ys = frexpl(y, &ey);
@@
-194,31
+195,41
@@
long double fmal(long double x, long double y, long double z)
* modes other than FE_TONEAREST are painful.
*/
if (spread < -LDBL_MANT_DIG) {
* modes other than FE_TONEAREST are painful.
*/
if (spread < -LDBL_MANT_DIG) {
+#ifdef FE_INEXACT
feraiseexcept(FE_INEXACT);
feraiseexcept(FE_INEXACT);
+#endif
+#ifdef FE_UNDERFLOW
if (!isnormal(z))
feraiseexcept(FE_UNDERFLOW);
if (!isnormal(z))
feraiseexcept(FE_UNDERFLOW);
+#endif
switch (oround) {
switch (oround) {
- case FE_TONEAREST:
+ default: /* FE_TONEAREST */
return (z);
return (z);
+#ifdef FE_TOWARDZERO
case FE_TOWARDZERO:
if (x > 0.0 ^ y < 0.0 ^ z < 0.0)
return (z);
else
return (nextafterl(z, 0));
case FE_TOWARDZERO:
if (x > 0.0 ^ y < 0.0 ^ z < 0.0)
return (z);
else
return (nextafterl(z, 0));
+#endif
+#ifdef FE_DOWNWARD
case FE_DOWNWARD:
if (x > 0.0 ^ y < 0.0)
return (z);
else
return (nextafterl(z, -INFINITY));
case FE_DOWNWARD:
if (x > 0.0 ^ y < 0.0)
return (z);
else
return (nextafterl(z, -INFINITY));
- default: /* FE_UPWARD */
+#endif
+#ifdef FE_UPWARD
+ case FE_UPWARD:
if (x > 0.0 ^ y < 0.0)
return (nextafterl(z, INFINITY));
else
return (z);
if (x > 0.0 ^ y < 0.0)
return (nextafterl(z, INFINITY));
else
return (z);
+#endif
}
}
if (spread <= LDBL_MANT_DIG * 2)
}
}
if (spread <= LDBL_MANT_DIG * 2)
- zs =
ldexp
l(zs, -spread);
+ zs =
scalbn
l(zs, -spread);
else
zs = copysignl(LDBL_MIN, zs);
else
zs = copysignl(LDBL_MIN, zs);
@@
-244,23
+255,25
@@
long double fmal(long double x, long double y, long double z)
*/
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) {
/*
* There is no need to worry about double rounding in directed
* rounding modes.
}
if (oround != FE_TONEAREST) {
/*
* There is no need to worry about double rounding in directed
* rounding modes.
+ * TODO: underflow is not raised correctly, example in downward rounding:
+ * fmal(0x1.0000000001p-16000L, 0x1.0000000001p-400L, -0x1p-16440L)
*/
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
);
else
else
- return
(add_and_denormalize(r.hi, adj, spread)
);
+ return
add_and_denormalize(r.hi, adj, spread
);
}
#endif
}
#endif