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.

Firearm

In compilers that support You can force the compiler to adhere to IEEE rules by using the appropriate intrinsic module. However, gfortran does not support that. You have to live with the fact that fast-math can be unsafe in this regard. – Hiawatha

When you turn on fast-math, you are promising the compiler that it can freely pretend that NaNs do not exist. This means that computations that would ordinarily produce NaN may not, and that the compiler can optimize away any code that checks for NaNs (since you promised that they don’t exist!), which makes isnan essentially useless. Your question is the equivalent of speeding on a curvy mountain road at night wearing a blindfold, and worrying about your busted tail light. – Bree

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..

Ivett

Mh have you checked that it actually always works with numbers that come out of modules that are compiled with -ffast-math? If so could you explain why? I thought the important information is in the number itself not in the checking routine. – Firearm

@Firearm No, the point is in the checking routine. The -ffast-math simply optimizes the check to false. It is very simple. – Hiawatha

Answer 2

In case someone is still looking for an answer, I find that the following works with gfortran -O3 -Ofast -ffast-math (gfortran version 9.3.0):

is_nan = (.not. (x <= huge(x) .and. x >= -huge(x))) .and. (.not. abs(x) > huge(x))

Indeed, I have posted on GitHub some functions for checking Inf/NaN, which aim to work even when compilers are invoked with aggressive optimization flags, such as gfortran -Ofast. Some simple tests are also included in the GitHub repository.

The ieee_is_nan included in ieee_arithmetic of gfortran 9.3.0 does not work with aggressive optimization flags like -Ofast (although it is not surprising). In addition, some compilers (gfortran 9.3.0, ifort 21.0, and nagfor 7.0) are buggy concerning the return type of ieee_is_nan when some special compilation options are imposed, as has been discussed here on Fortran Discourse.

For your convenience, I copy-paste the main code below. More details can be found on GitHub.

! consts.f90
module consts_mod

implicit none
private
public :: SP, DP

integer, parameter :: SP = kind(0.0)
integer, parameter :: DP = kind(0.0D0)

end module consts_mod

! infnan.f90
module infnan_mod

use consts_mod, only : SP, DP
implicit none
private
public :: is_nan, is_finite, is_inf, is_posinf, is_neginf


interface is_nan
    module procedure is_nan_sp, is_nan_dp
end interface is_nan
 
interface is_finite
    module procedure is_finite_sp, is_finite_dp
end interface is_finite

interface is_posinf
    module procedure is_posinf_sp, is_posinf_dp
end interface is_posinf

interface is_neginf
    module procedure is_neginf_sp, is_neginf_dp
end interface is_neginf

interface is_inf
    module procedure is_inf_sp, is_inf_dp
end interface is_inf


contains


elemental pure function is_nan_sp(x) result(y)
implicit none
real(SP), intent(in) :: x
logical :: y
y = (.not. (x <= huge(x) .and. x >= -huge(x))) .and. (.not. abs(x) > huge(x))
end function is_nan_sp

elemental pure function is_nan_dp(x) result(y)
implicit none
real(DP), intent(in) :: x
logical :: y
y = (.not. (x <= huge(x) .and. x >= -huge(x))) .and. (.not. abs(x) > huge(x))
end function is_nan_dp


elemental pure function is_finite_sp(x) result(y)
implicit none
real(SP), intent(in) :: x
logical :: y
y = (x <= huge(x) .and. x >= -huge(x))
end function is_finite_sp

elemental pure function is_finite_dp(x) result(y)
implicit none
real(DP), intent(in) :: x
logical :: y
y = (x <= huge(x) .and. x >= -huge(x))
end function is_finite_dp


elemental pure function is_inf_sp(x) result(y)
implicit none
real(SP), intent(in) :: x
logical :: y
y = (abs(x) > huge(x))
end function is_inf_sp

elemental pure function is_inf_dp(x) result(y)
implicit none
real(DP), intent(in) :: x
logical :: y
y = (abs(x) > huge(x))
end function is_inf_dp


elemental pure function is_posinf_sp(x) result(y)
implicit none
real(SP), intent(in) :: x
logical :: y
y = (abs(x) > huge(x)) .and. (x > 0)
end function is_posinf_sp

elemental pure function is_posinf_dp(x) result(y)
implicit none
real(DP), intent(in) :: x
logical :: y
y = (abs(x) > huge(x)) .and. (x > 0)
end function is_posinf_dp


elemental pure function is_neginf_sp(x) result(y)
implicit none
real(SP), intent(in) :: x
logical :: y
y = (abs(x) > huge(x)) .and. (x < 0)
end function is_neginf_sp

elemental pure function is_neginf_dp(x) result(y)
implicit none
real(DP), intent(in) :: x
logical :: y
y = (abs(x) > huge(x)) .and. (x < 0)
end function is_neginf_dp


end module infnan_mod

Strobe

"The ieee_is_nan included in ieee_arithmetic of gfortran 9.3.0 does not work with aggressive optimization flags like -Ofast." That isn't really surprising given the explanation in my answer. The check is, most likely after inlining, optimized away in the very same way as a manual x/=x or an equivalent. I wonder why my answer had to be downvoted (not accusing you personally). – Hiawatha

Functions that are working regardless of optimization flags are of course very useful. They just require more arithmetic operations and the performance has to be tested. The approach in my answer prevents inlining and might be even slower. – Hiawatha

I don't understand, isn't it better from the point of view of performance or reliability, isn't it better to simply compare the argument with constants (or preferably initialized variables)? For example: y = (x .eq. SP_signaling_NaN) .or. (x .eq. SP_quiet_NaN). – Sorbose

We can not check whether x is NaN by testing whether the equality x == NaN holds --- it does not, which is because of the very definition of NaN. – Strobe

Answer 3

I was facing the same problem in gfortran-4.7 and started to experiment with some expressions. For my Test Cases this one returned true if x is a NaN:

check = ((x*2d0==x).AND.(
         (x<-epsilon(1.d0)).OR.
         (x>+epsilon(1.d0))))

I've checked only double precision values and -O2 -ffast-math, -O3 -ffast-math options. Note that this returns .false. if using no optimization flags, thus you would have to combine it with

isnan(x) .OR. check

Extramarital

That is an interesting hack but not really a solution if you want to support multiple compilers and versions and precision types :). – Firearm