Fix NaN detection with -Ofast

* isnan and isfinite does not work when building with -Ofast
* The issue has always been there but it's now detected by clang >= 18
* This version of isfinite may be slower but we only use it for debug

src/compression.cpp

@@ -44,18 +44,6 @@
 #include "cluster_assembly_function.hpp"
 #include "random_pivot_manager.hpp"
 
-#ifdef _MSC_VER
-// Intel compiler defines isnan in global namespace
-// MSVC defines _isnan
-# ifndef __INTEL_COMPILER
-#  define isnan _isnan
-# endif
-#elif __GLIBC__ == 2 && __GLIBC_MINOR__ < 23
-// https://sourceware.org/bugzilla/show_bug.cgi?id=19439
-#elif __cplusplus >= 201103L || !defined(__GLIBC__)
-using std::isnan;
-#endif
-
 using std::vector;
 using std::min;
 using std::max;
@@ -764,7 +752,7 @@ template<typename T> RkMatrix<typename Types<T>::dp>* compressOneStratum(
 
     // If I meet a NaN, I save & leave
     // TODO : improve this behaviour
-    if (isnan(approxNorm)) {
+    if (!swIsFinite(approxNorm)) {
       rkFull->toFile("Rk");
       full->toFile("Full");
       HMAT_ASSERT(false);

src/data_types.cpp

@@ -21,6 +21,8 @@
 */
 
 #include "data_types.hpp"
+#include <cstdint>
+#include <cstring>
 
 namespace hmat {
 
@@ -38,4 +40,19 @@ template<> const int Constants<D_t>::code = 1;
 template<> const int Constants<C_t>::code = 2;
 template<> const int Constants<Z_t>::code = 3;
 
+template<> bool swIsFinite(double x) {
+  uint64_t bits;
+  // casting break strict aliasing rules and std::bit_cast is C++20 only
+  std::memcpy(&bits, &x, sizeof(double));
+  // FIXME: This trick fool gcc up to 14.x and clang up to 18.x. Clang 19 will
+  // optimize and always return true when using -ffast-math
+  return (bits & 0x7ff0000000000000) != 0x7ff0000000000000;
+}
+
+template<> bool swIsFinite(float x) {
+  uint32_t bits;
+  std::memcpy(&bits, &x, sizeof(float));
+  return (bits & 0x7f800000) != 0x7f800000;
+}
+
 }  // end namespace hmat

src/data_types.hpp

@@ -146,6 +146,11 @@ inline Z_t conj(const Z_t x) {
   return std::conj(x);
 }
 
+/** Same as std::isfinite but does not rely on FPU so work with -Ofast */
+template<typename T> bool swIsFinite(T);
+template<> bool swIsFinite(double x);
+template<> bool swIsFinite(float x);
+
 }  // end namespace hmat
 
 #endif

src/scalar_array.cpp

@@ -612,7 +612,7 @@ template<typename T, typename std::enable_if<hmat::Types<T>::IS_REAL::value, T*>
 void checkNanSFINAE(const ScalarArray<T>* m) {
   for (int col = 0; col < m->cols; col++) {
     for (int row = 0; row < m->rows; row++) {
-      HMAT_ASSERT(std::isfinite(m->get(row, col)));
+      HMAT_ASSERT(swIsFinite(m->get(row, col)));
     }
   }
 }
@@ -621,8 +621,8 @@ template<typename T, typename std::enable_if<!hmat::Types<T>::IS_REAL::value, T*
 void checkNanSFINAE(const ScalarArray<T>* m) {
   for (int col = 0; col < m->cols; col++) {
     for (int row = 0; row < m->rows; row++) {
-      HMAT_ASSERT(std::isfinite(m->get(row, col).real()));
-      HMAT_ASSERT(std::isfinite(m->get(row, col).imag()));
+      HMAT_ASSERT(swIsFinite(m->get(row, col).real()));
+      HMAT_ASSERT(swIsFinite(m->get(row, col).imag()));
     }
   }
 }