math: add fma TODO comments about the underflow issue

The underflow exception is not raised correctly in some
cornercases (see previous fma commit), added comments
with examples for fmaf, fmal and non-x86 fma.

In fmaf store the result before returning so it has the
correct precision when FLT_EVAL_METHOD!=0

diff --git a/src/math/fma.c b/src/math/fma.c
index 89def795943573d1568c5328ee95b2c9efc27b36..850a4a6c7f0dcc164521082be8fef0fac986c11f 100644 (file)
--- a/src/math/fma.c
+++ b/src/math/fma.c
@@ -441,6 +441,8 @@ double fma(double x, double y, double z)
                /*
                 * 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;

diff --git a/src/math/fmaf.c b/src/math/fmaf.c
index a1c7f4f83df1e3bb0f6c67c37765345db5be4feb..745ee393c553d171d85b4cfcce06eccfcb4378df 100644 (file)
--- a/src/math/fmaf.c
+++ b/src/math/fmaf.c
@@ -49,7 +49,14 @@ float fmaf(float x, float y, float z)
                (hr & 0x7ff00000) == 0x7ff00000 ||  /* NaN */
                result - xy == z ||                 /* exact */
                fegetround() != FE_TONEAREST)       /* not round-to-nearest */
-               return (result);
+       {
+               /*
+               TODO: underflow is not raised correctly, example in
+               downward rouding: fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)
+               */
+               z = result;
+               return z;
+       }
 
        /*
         * If result is inexact, and exactly halfway between two float values,
@@ -63,5 +70,6 @@ float fmaf(float x, float y, float z)
        fesetround(FE_TONEAREST);
        if (result == adjusted_result)
                SET_LOW_WORD(adjusted_result, lr + 1);
-       return (adjusted_result);
+       z = adjusted_result;
+       return z;
 }

diff --git a/src/math/fmal.c b/src/math/fmal.c
index ccbe434d0eae6b883fc07e48bc3b5b79f9c971ff..87e30fcf65394f71b34027081feb9cab0f013ca4 100644 (file)
--- a/src/math/fmal.c
+++ b/src/math/fmal.c
@@ -262,6 +262,8 @@ long double fmal(long double x, long double y, long double z)
                /*
                 * 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;