Provide always-working isnan etc.

Certain math optimization options (e.g. -ffast-math) tell the compiler
that IEEE floating point numbers such as inf and nan do not need to be
handled correctly (in the sense specified by the IEEE standard). This
greatly improves floating-point speed and is commonly used in
numerical HPC applications.

However, since compilers then don't need to handle inf and nan
correctly, they have begun to optimise isnan(x) to simply returning
false all the time. This improves speed (since the check does not
actually need to occur) and reduces code size (since the nan-handling
if branches can be omitted). Of course, this makes it then impossible
to actually check for nan by calling isnan.

Currently, e.g. g++ performs this optimisation, whereas gcc does not.
Things vary with other compilers. In the future, with link-time
optimisations, I expect other compilers to follow g++.

This patch provides functions CCTK_IEEE_isnan etc. that always check
for nan, independent of the chosen optimisation flags.


git-svn-id: http://svn.cactuscode.org/flesh/trunk@5107 17b73243-c579-4c4c-a9d2-2d5706c11dac

2 files changed, 91 insertions, 0 deletions

diff --git a/src/include/cctk_Math.h b/src/include/cctk_Math.h
index e51c991f..5bc29b59 100644
--- a/src/include/cctk_Math.h
+++ b/src/include/cctk_Math.h
@@ -41,6 +41,13 @@ int CCTK_isnan(double x);
 int CCTK_isnormal(double x);
 int CCTK_signbit(double x);
 
+  /* int CCTK_IEEE_fpclassify(double x); */
+int CCTK_IEEE_isfinite(double x);
+int CCTK_IEEE_isinf(double x);
+int CCTK_IEEE_isnan(double x);
+int CCTK_IEEE_isnormal(double x);
+int CCTK_IEEE_signbit(double x);
+
 #ifdef __cplusplus
 }   
 #endif
diff --git a/src/util/Math.c b/src/util/Math.c
index 7ce8d680..dbc146ad 100644
--- a/src/util/Math.c
+++ b/src/util/Math.c
@@ -10,6 +10,8 @@
  @@*/
 
 #include <math.h>
+#include <stdint.h>
+#include <string.h>
 #include <cctk_Config.h>
 #include <cctk_Math.h>
 
@@ -131,3 +133,85 @@ int CCTK_FNAME(CCTK_signbit)(double *x)
 {
   return CCTK_signbit(*x);
 }
+
+
+
+int CCTK_IEEE_isfinite(double x)
+{
+  uint64_t u;
+  memcpy(&u, &x, sizeof u);
+  uint64_t e = (u & UINT64_C(0x7ff0000000000000)) >> 52;
+  return e != UINT64_C(0x7ff);
+}
+
+int CCTK_FNAME(CCTK_IEEE_isfinite)(double *x);
+int CCTK_FNAME(CCTK_IEEE_isfinite)(double *x)
+{
+  return CCTK_IEEE_isfinite(*x);
+}
+
+
+
+int CCTK_IEEE_isinf(double x)
+{
+  uint64_t u;
+  memcpy(&u, &x, sizeof u);
+  uint64_t e = (u & UINT64_C(0x7ff0000000000000)) >> 52;
+  uint64_t m = u & UINT64_C(0x000fffffffffffff);
+  return e == UINT64_C(0x7ff) && m == 0;
+}
+
+int CCTK_FNAME(CCTK_IEEE_isinf)(double *x);
+int CCTK_FNAME(CCTK_IEEE_isinf)(double *x)
+{
+  return CCTK_IEEE_isinf(*x);
+}
+
+
+
+int CCTK_IEEE_isnan(double x)
+{
+  uint64_t u;
+  memcpy(&u, &x, sizeof u);
+  uint64_t e = (u & UINT64_C(0x7ff0000000000000)) >> 52;
+  uint64_t m = u & UINT64_C(0x000fffffffffffff);
+  return e == UINT64_C(0x7ff) && m != 0;
+}
+
+int CCTK_FNAME(CCTK_IEEE_isnan)(double *x);
+int CCTK_FNAME(CCTK_IEEE_isnan)(double *x)
+{
+  return CCTK_IEEE_isnan(*x);
+}
+
+
+
+int CCTK_IEEE_isnormal(double x)
+{
+  uint64_t u;
+  memcpy(&u, &x, sizeof u);
+  uint64_t e = (u & UINT64_C(0x7ff0000000000000)) >> 52;
+  return e != 0 && e != UINT64_C(0x7ff);
+}
+
+int CCTK_FNAME(CCTK_IEEE_isnormal)(double *x);
+int CCTK_FNAME(CCTK_IEEE_isnormal)(double *x)
+{
+  return CCTK_IEEE_isnormal(*x);
+}
+
+
+
+int CCTK_IEEE_signbit(double x)
+{
+  uint64_t u;
+  memcpy(&u, &x, sizeof u);
+  uint64_t s = (u & UINT64_C(0x8000000000000000)) >> 63;
+  return s;
+}
+
+int CCTK_FNAME(CCTK_IEEE_signbit)(double *x);
+int CCTK_FNAME(CCTK_IEEE_signbit)(double *x)
+{
+  return CCTK_IEEE_signbit(*x);
+}