| Index: fusl/src/math/fmaf.c
|
| diff --git a/fusl/src/math/fmaf.c b/fusl/src/math/fmaf.c
|
| index 4ebe855204d649f76be1445ad2e27f7e97cd8e5b..b3ecd37bf8b08d140fc705feb0959dbce55d5a47 100644
|
| --- a/fusl/src/math/fmaf.c
|
| +++ b/fusl/src/math/fmaf.c
|
| @@ -37,58 +37,60 @@
|
| * direct double-precision arithmetic suffices, except where double
|
| * rounding occurs.
|
| */
|
| -float fmaf(float x, float y, float z)
|
| -{
|
| - PRAGMA_STDC_FENV_ACCESS_ON
|
| - double xy, result;
|
| - union {double f; uint64_t i;} u;
|
| - int e;
|
| +float fmaf(float x, float y, float z) {
|
| + PRAGMA_STDC_FENV_ACCESS_ON
|
| + double xy, result;
|
| + union {
|
| + double f;
|
| + uint64_t i;
|
| + } u;
|
| + int e;
|
|
|
| - xy = (double)x * y;
|
| - result = xy + z;
|
| - u.f = result;
|
| - e = u.i>>52 & 0x7ff;
|
| - /* Common case: The double precision result is fine. */
|
| - if ((u.i & 0x1fffffff) != 0x10000000 || /* not a halfway case */
|
| - e == 0x7ff || /* NaN */
|
| - result - xy == z || /* exact */
|
| - fegetround() != FE_TONEAREST) /* not round-to-nearest */
|
| - {
|
| - /*
|
| - underflow may not be raised correctly, example:
|
| - fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)
|
| - */
|
| + xy = (double)x * y;
|
| + result = xy + z;
|
| + u.f = result;
|
| + e = u.i >> 52 & 0x7ff;
|
| + /* Common case: The double precision result is fine. */
|
| + if ((u.i & 0x1fffffff) != 0x10000000 || /* not a halfway case */
|
| + e == 0x7ff || /* NaN */
|
| + result - xy == z || /* exact */
|
| + fegetround() != FE_TONEAREST) /* not round-to-nearest */
|
| + {
|
| +/*
|
| +underflow may not be raised correctly, example:
|
| +fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)
|
| +*/
|
| #if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
|
| - if (e < 0x3ff-126 && e >= 0x3ff-149 && fetestexcept(FE_INEXACT)) {
|
| - feclearexcept(FE_INEXACT);
|
| - /* TODO: gcc and clang bug workaround */
|
| - volatile float vz = z;
|
| - result = xy + vz;
|
| - if (fetestexcept(FE_INEXACT))
|
| - feraiseexcept(FE_UNDERFLOW);
|
| - else
|
| - feraiseexcept(FE_INEXACT);
|
| - }
|
| + if (e < 0x3ff - 126 && e >= 0x3ff - 149 && fetestexcept(FE_INEXACT)) {
|
| + feclearexcept(FE_INEXACT);
|
| + /* TODO: gcc and clang bug workaround */
|
| + volatile float vz = z;
|
| + result = xy + vz;
|
| + if (fetestexcept(FE_INEXACT))
|
| + feraiseexcept(FE_UNDERFLOW);
|
| + else
|
| + feraiseexcept(FE_INEXACT);
|
| + }
|
| #endif
|
| - z = result;
|
| - return z;
|
| - }
|
| + z = result;
|
| + return z;
|
| + }
|
|
|
| - /*
|
| - * If result is inexact, and exactly halfway between two float values,
|
| - * we need to adjust the low-order bit in the direction of the error.
|
| - */
|
| +/*
|
| + * If result is inexact, and exactly halfway between two float values,
|
| + * we need to adjust the low-order bit in the direction of the error.
|
| + */
|
| #ifdef FE_TOWARDZERO
|
| - fesetround(FE_TOWARDZERO);
|
| + fesetround(FE_TOWARDZERO);
|
| #endif
|
| - volatile double vxy = xy; /* XXX work around gcc CSE bug */
|
| - double adjusted_result = vxy + z;
|
| - fesetround(FE_TONEAREST);
|
| - if (result == adjusted_result) {
|
| - u.f = adjusted_result;
|
| - u.i++;
|
| - adjusted_result = u.f;
|
| - }
|
| - z = adjusted_result;
|
| - return z;
|
| + volatile double vxy = xy; /* XXX work around gcc CSE bug */
|
| + double adjusted_result = vxy + z;
|
| + fesetround(FE_TONEAREST);
|
| + if (result == adjusted_result) {
|
| + u.f = adjusted_result;
|
| + u.i++;
|
| + adjusted_result = u.f;
|
| + }
|
| + z = adjusted_result;
|
| + return z;
|
| }
|
|
|