Is it possible to make isnan() work in gfortran -O3 -ffast-math?

Question

I would like to compile a program with gfortran and -O3 -ffast-math enabled, since it gives a nice performance boost. I was rather confused, that gfortran's isnan() catched some NaN's but not all of them. After reading

Checking if a double (or float) is NaN in C++
how do I make a portable isnan/isinf function
Negative NaN is not a NaN?

I am under the impression that people are able to check for NaN's in C via bit-fiddling even with fast-math enabled. However, this puzzles me since fast-math

can result in incorrect output for programs that depend on an exact implementation of IEEE or ISO rules/specifications for math functions.

According to the man page of gcc 4.7.2. So how do you know which bit to check, if the numbers are not represented according to IEEE standard? And if you know it, how would you implement it in Fortran 95/03/08?

Don't bother to post (x \= x) or simlar solutions which depend on IEEE rules. They give the same result as isnan(). I am also aware of -ffpe-trap=invalid,zero,overflow, but don't want to stop the program. If it helps, my OS is 64-bit LinuxMint 14. If it is not possible in Fortran, a waterproof C solution would also be nice.

Answer 1

First I would point out that gfortran 4.9 supports the IEEE_arithmetic module. However, the checks in the procedures in that module can be optimized away and the same logic applies.

However, I could not depend GCC 4.9in 2014, it was too fresh. I used the following.

When I had to use x/=x I moved the check x/=x to a procedure which is compiled without -ffast-math and without link time optimizations:

module my_is_nan_mod
  !poor man's replacement for the intrinsic module
  !do not use if the compiler supports the F2003 version
  !make sure not to use fast math optimizations when compiling
  use iso_fortran_env
  
  !the common extension isnan() may actually fail with fast math optimizations

  
  interface my_is_nan
    module procedure my_is_nan_real32
    module procedure my_is_nan_real64
  end interface

contains
  logical function my_is_nan_real32(x) result(res)
    real(real32), intent(in) :: x

    res = x /= x
  end function

  logical elemental function my_is_nan_real64(x) result(res)
    real(real64), intent(in) :: x

    res = x /= x
  end function

end module

It is in a separate file, which is then compiled without -Ofast, --ffast-math and without -flto. Beware the lack of inlining can cause serious performance decrease.

With -ffast-math the compiler sees x/=x, decides that it cannot ever be true and optimizes the expression to just .false..