nextafterf: trying to correct FP exception handling

diff --git a/docs/probe_math_exception.c b/docs/probe_math_exception.c
new file mode 100644 (file)
index 0000000..dbeccc5
--- /dev/null
+++ b/docs/probe_math_exception.c
@@ -0,0 +1,41 @@
+/* Small test program for probing how various math functions
+ * with specific operands set floating point exceptions
+ */
+
+#define _ISOC99_SOURCE 1
+#define _GNU_SOURCE    1
+
+#include <math.h>
+#include <fenv.h>
+#include <stdio.h>
+
+int main(int argc, char **argv)
+{
+       float infF = HUGE_VALF * 2;
+
+       feclearexcept(FE_ALL_EXCEPT);
+
+//     printf("%.40e\n", 1.0 / 0.0); // FE_DIVBYZERO
+//     printf("%.40e\n", nextafterf(HUGE_VALF, infF)); // no exceptions in glibc 2.4
+
+#define PREX(ex) do { if (fetestexcept(ex)) printf(#ex); } while(0)
+#ifdef FE_INEXACT
+       PREX(FE_INEXACT);
+#endif
+#ifdef FE_DIVBYZERO
+       PREX(FE_DIVBYZERO);
+#endif
+#ifdef FE_UNDERFLOW
+       PREX(FE_UNDERFLOW);
+#endif
+#ifdef FE_OVERFLOW
+       PREX(FE_OVERFLOW);
+#endif
+#ifdef FE_INVALID
+       PREX(FE_INVALID);
+#endif
+       if (fetestexcept(FE_ALL_EXCEPT))
+               printf("\n");
+       printf("done\n");
+       return 0;
+}

diff --git a/libm/math_private.h b/libm/math_private.h
index bdd0aba..0601f2c 100644 (file)
--- a/libm/math_private.h
+++ b/libm/math_private.h
@@ -187,4 +187,59 @@ extern double __kernel_cos (double,double) attribute_hidden;
 extern double __kernel_tan (double,double,int) attribute_hidden;
 extern int    __kernel_rem_pio2 (double*,double*,int,int,int,const int*) attribute_hidden;
 
+/*
+ * math_opt_barrier(x): force expression x to be evaluated and put into
+ * a floating point register or memory. This macro returns the value.
+ *
+ * math_force_eval(x): force expression x to be evaluated and put into
+ * a floating point register or memory *of the appropriate size*.
+ * This forces floating point flags to be set correctly
+ * (for example, when float value is overflowing, but FPU registers
+ * are wide enough to "hide" this).
+ */
+#if defined(__i386__)
+#define math_opt_barrier(x) ({ \
+       __typeof(x) __x = (x); \
+       /* "t": load x into top-of-stack fpreg */ \
+       __asm ("" : "=t" (__x) : "0" (__x)); \
+       __x; \
+})
+#define math_force_eval(x) do {        \
+       if (sizeof(x) <= sizeof(double)) \
+               /* "m": store x into a memory location */ \
+               __asm __volatile ("" : : "m" (x)); \
+       else /* long double */ \
+               /* "f": load x into (any) fpreg */ \
+               __asm __volatile ("" : : "f" (x)); \
+} while (0)
+#endif
+
+#if defined(__x86_64__)
+#define math_opt_barrier(x) ({ \
+       __typeof(x) __x = (x); \
+       if (sizeof(x) <= sizeof(double)) \
+               /* "x": load into XMM SSE register */ \
+               __asm ("" : "=x" (__x) : "0" (__x)); \
+       else /* long double */ \
+               /* "t": load x into top-of-stack fpreg */ \
+               __asm ("" : "=t" (__x) : "0" (__x)); \
+       __x; \
+})
+#define math_force_eval(x) do { \
+       if (sizeof(x) <= sizeof(double)) \
+               /* "x": load into XMM SSE register */ \
+               __asm __volatile ("" : : "x" (x)); \
+       else /* long double */ \
+               /* "f": load x into (any) fpreg */ \
+               __asm __volatile ("" : : "f" (x)); \
+} while (0)
+#endif
+
+/* Default implementation forces store to a memory location */
+#ifndef math_opt_barrier
+#define math_opt_barrier(x) ({ __typeof(x) __x = (x); __asm ("" : "+m" (__x)); __x; })
+#define math_force_eval(x)  __asm __volatile ("" : : "m" (x))
+#endif
+
+
 #endif /* _MATH_PRIVATE_H_ */

diff --git a/libm/s_nextafterf.c b/libm/s_nextafterf.c
index 8dee00f..5fc44e3 100644 (file)
--- a/libm/s_nextafterf.c
+++ b/libm/s_nextafterf.c
@@ -16,11 +16,6 @@
 #include "math.h"
 #include "math_private.h"
 
-#ifndef math_opt_barrier
-# define math_opt_barrier(x) ({ __typeof (x) __x = x; __asm ("" : "+m" (__x)); __x; })
-# define math_force_eval(x)  __asm __volatile ("" : : "m" (x))
-#endif
-
 float nextafterf(float x, float y)
 {
        int32_t hx, hy, ix, iy;
@@ -38,12 +33,12 @@ float nextafterf(float x, float y)
                return y;
 
        if (ix == 0) { /* x == 0? */
-               float u;
+// glibc 2.4 does not seem to set underflow?
+//             float u;
                /* return +-minsubnormal */
                SET_FLOAT_WORD(x, (hy & 0x80000000) | 1);
-               u = math_opt_barrier(x);
-               u = u * u;
-               math_force_eval(u); /* raise underflow flag */
+//             u = x * x; /* raise underflow flag */
+//             math_force_eval(u);
                return x;
        }
 
@@ -63,8 +58,6 @@ float nextafterf(float x, float y)
        hy = hx & 0x7f800000;
        if (hy >= 0x7f800000) {
                x = x + x; /* overflow */
-//??           if (FLT_EVAL_METHOD != 0)
-//                     asm ("" : "+m"(x));
                return x; /* overflow */
        }
        if (hy < 0x00800000) {