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; |
} |