*/
#include <fenv.h>
-#include "libm.h"
+#include <math.h>
+#include <stdint.h>
/*
* Fused multiply-add: Compute x * y + z with a single rounding error.
*/
float fmaf(float x, float y, float z)
{
+ #pragma STDC FENV_ACCESS ON
double xy, result;
- uint32_t hr, lr;
+ union {double f; uint64_t i;} u;
+ int e;
xy = (double)x * y;
result = xy + z;
- EXTRACT_WORDS(hr, lr, result);
+ u.f = result;
+ e = u.i>>52 & 0x7ff;
/* Common case: The double precision result is fine. */
- if ((lr & 0x1fffffff) != 0x10000000 || /* not a halfway case */
- (hr & 0x7ff00000) == 0x7ff00000 || /* NaN */
- result - xy == z || /* exact */
+ if ((u.i & 0x1fffffff) != 0x10000000 || /* not a halfway case */
+ e == 0x7ff || /* NaN */
+ (result - xy == z && result - z == xy) || /* exact */
fegetround() != FE_TONEAREST) /* not round-to-nearest */
- return (result);
+ {
+ /*
+ 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);
+ }
+#endif
+ 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.
*/
-#ifdef 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)
- SET_LOW_WORD(adjusted_result, lr + 1);
- return (adjusted_result);
+ double err;
+ int neg = u.i >> 63;
+ if (neg == (z > xy))
+ err = xy - result + z;
+ else
+ err = z - result + xy;
+ if (neg == (err < 0))
+ u.i++;
+ else
+ u.i--;
+ z = u.f;
+ return z;
}