/*
assume (long double)(hi+lo) == hi
-return an adjusted hi so that rounding it to double is correct
+return an adjusted hi so that rounding it to double (or less) precision is correct
*/
static long double adjust(long double hi, long double lo)
{
ulo.x = lo;
if (uhi.bits.s == ulo.bits.s)
uhi.bits.m++;
- else
+ else {
uhi.bits.m--;
+ /* handle underflow and take care of ld80 implicit msb */
+ if (uhi.bits.m == (uint64_t)-1/2) {
+ uhi.bits.m *= 2;
+ uhi.bits.e--;
+ }
+ }
return uhi.x;
}
+/* adjusted add so the result is correct when rounded to double (or less) precision */
static long double dadd(long double x, long double y)
{
add(&x, &y, x, y);
return adjust(x, y);
}
+/* adjusted mul so the result is correct when rounded to double (or less) precision */
static long double dmul(long double x, long double y)
{
mul(&x, &y, x, y);
double fma(double x, double y, double z)
{
+ #pragma STDC FENV_ACCESS ON
long double hi, lo1, lo2, xy;
int round, ez, exy;
} else if (ez > exy - 12) {
add(&hi, &lo2, xy, z);
if (hi == 0) {
+ /*
+ xy + z is 0, but it should be calculated with the
+ original rounding mode so the sign is correct, if the
+ compiler does not support FENV_ACCESS ON it does not
+ know about the changed rounding mode and eliminates
+ the xy + z below without the volatile memory access
+ */
+ volatile double z_;
fesetround(round);
- /* make sure that the sign of 0 is correct */
- return (xy + z) + lo1;
+ z_ = z;
+ return (xy + z_) + lo1;
}
} else {
/*
hi = xy;
lo2 = z;
}
+ /*
+ the result is stored before return for correct precision and exceptions
+
+ one corner case is when the underflow flag should be raised because
+ the precise result is an inexact subnormal double, but the calculated
+ long double result is an exact subnormal double
+ (so rounding to double does not raise exceptions)
+
+ in nearest rounding mode dadd takes care of this: the last bit of the
+ result is adjusted so rounding sees an inexact value when it should
+
+ in non-nearest rounding mode fenv is used for the workaround
+ */
fesetround(round);
if (round == FE_TONEAREST)
- return dadd(hi, dadd(lo1, lo2));
- return hi + (lo1 + lo2);
+ z = dadd(hi, dadd(lo1, lo2));
+ else {
+#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
+ int e = fetestexcept(FE_INEXACT);
+ feclearexcept(FE_INEXACT);
+#endif
+ z = hi + (lo1 + lo2);
+#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
+ if (getexp(z) < 0x3fff-1022 && fetestexcept(FE_INEXACT))
+ feraiseexcept(FE_UNDERFLOW);
+ else if (e)
+ feraiseexcept(FE_INEXACT);
+#endif
+ }
+ return z;
}
#else
/* origin: FreeBSD /usr/src/lib/msun/src/s_fma.c */
*/
double fma(double x, double y, double z)
{
+ #pragma STDC FENV_ACCESS ON
double xs, ys, zs, adj;
struct dd xy, r;
int oround;
/*
* There is no need to worry about double rounding in directed
* rounding modes.
+ * TODO: underflow is not raised properly, example in downward rounding:
+ * fma(0x1.000000001p-1000, 0x1.000000001p-30, -0x1p-1066)
*/
fesetround(oround);
adj = r.lo + xy.lo;