Challenge: Testing Inf and NaN with `gfortran-13 -Ofast`

[Update 20230911: There is a potential solution suggested on GitHub, but it contains non-standard code. Thank you for giving your comments.]

Here is a challenge for those who are interested. It is a true challenge for me, but maybe I have overlooked something and it turns out trivial for you.

Could anyone write five functions is_nan, is_inf, is_posinf, is_neginf, and is_finite to test whether each entry of x and y in the following code is NaN, Infinity (+ or -), +Infinity, -Infinity, or finite? The functions must work correctly when the WHOLE program is compiled with

gfortran-13 -Ofast

Note that ieee_is_nan etc. from ieee_arithmetic do not work under this situation.

Code:

PINF_SP = ieee_value(1.0, ieee_positive_inf)
NINF_SP = ieee_value(1.0, ieee_negative_inf)
QNAN_SP = ieee_value(1.0, ieee_quiet_nan)
SNAN_SP = ieee_value(1.0, ieee_signaling_nan)
PNORMAL_SP = ieee_value(1.0, ieee_positive_normal)
NNORMAL_SP = ieee_value(1.0, ieee_negative_normal)
PDENORMAL_SP = ieee_value(1.0, ieee_positive_denormal)
NDENORMAL_SP = ieee_value(1.0, ieee_negative_denormal)
PZERO_SP = ieee_value(1.0, ieee_positive_zero)
NZERO_SP = ieee_value(1.0, ieee_negative_zero)

PINF_DP = ieee_value(1.0D0, ieee_positive_inf)
NINF_DP = ieee_value(1.0D0, ieee_negative_inf)
QNAN_DP = ieee_value(1.0D0, ieee_quiet_nan)
SNAN_DP = ieee_value(1.0D0, ieee_signaling_nan)
PNORMAL_DP = ieee_value(1.0D0, ieee_positive_normal)
NNORMAL_DP = ieee_value(1.0D0, ieee_negative_normal)
PDENORMAL_DP = ieee_value(1.0D0, ieee_positive_denormal)
NDENORMAL_DP = ieee_value(1.0D0, ieee_negative_denormal)
PZERO_DP = ieee_value(1.0D0, ieee_positive_zero)
NZERO_DP = ieee_value(1.0D0, ieee_negative_zero)

x = [PINF_SP, NINF_SP, &
    & QNAN_SP, SNAN_SP, &
    & PNORMAL_SP, NNORMAL_SP, &
    & PDENORMAL_SP, NDENORMAL_SP, &
    & PZERO_SP, NZERO_SP, &
    & huge(0.0), -huge(0.0), &
    & epsilon(0.0), -epsilon(0.0), &
    & tiny(0.0), -tiny(0.0), &
    & 0.0, 1.0, -1.0]

y = [PINF_DP, NINF_DP, &
    & QNAN_DP, SNAN_DP, &
    & PNORMAL_DP, NNORMAL_DP, &
    & PDENORMAL_DP, NDENORMAL_DP, &
    & PZERO_DP, NZERO_DP, &
    & huge(0.0D0), -huge(0.0D0), &
    & epsilon(0.0D0), -epsilon(0.0D0), &
    & tiny(0.0D0), -tiny(0.0D0), &
    & 0.0D0, 1.0D0, -1.0D0]

Motivation: I have a project that provides such functions, which are supposed to work correctly even if compilers are invoked with aggressive optimization flags. They work without any problem with gfortran 9, 10, 11, and 12 (in addition to many other compilers), but fail with gfortran-13. If you are interested, you may clone my project and do

cd test && make gtest

which should fail if your gfortran is gfortran 13. You may also test ifort with make itest.

Below is the Fortran file that I use to conduct the above-mentioned test. Note that it depends on other files in the project.

Thanks.

N.B.: It is well known why naive testing of Inf and NaN can fail when compilers are invoked with aggressive optimization options. To be clear, the “why” is NOT the topic of this thread. The topic is “how”. Thanks.

3 Likes

You could try this one: Brad Richardson / exceptional_numbers · GitLab

Thank you @everythingfunctional .

Our implementations of these functions share similar ideas. See https://github.com/equipez/infnan/blob/main/infnan.f90 and https://github.com/equipez/infnan/blob/main/inf.f90.

I made a quick test using your code, but is_nan(QNAN_SP) and is_nan(SNAN_SP) both return .false. if the code is compiled with gfortran-13 -Ofast.

Thank you.

:frowning: that’s too bad.

Edit: to confirm, the library’s own test suite fails the following cases when compiled with -Ofast

Test that
    is_infinity
        is true for +inf
            Expected to be true
        is true for -inf
            Expected to be true
    is_nan
        is true for nan
            Expected to be true
    is_zero
        is false for number smaller than tiny
            Expected to not be true
        is false for nan
            Expected to not be true

This seems to be a very challenging problem.

Hi @zaikunzhang I tried tweaking the functions but to no avail with gfortran 13 -Ofast, but with -O3 everything works. Just wondering, do you get such big differences between -O3 and -Ofast? (Pure curiosity as there seems to be indeed a regression in gfortran, versions <13 do work as you said)

1 Like

Is it possibly useful to compile the source file containing the custom test functions like is_nan() with -O3 only (i.e., no -ffast-math or -Ofast), and then link the obtained object file with other object files compiles with -Ofast etc?

FWIW, in March this year, I had experienced a severe / very troublesome bug with the behavior of literal NaN and its use as a sentinel value in my codes, which was caused by the “unexpected” behavior of builtin isnan() or ieee_is_nan() in combination with an option involving fast-math. Since then, I completely stopped using both NaN and fast-math in my codes (because I cannot guarantee them to work in a portable way for different compilers/options/environments). I think the current situation of NaN in Fortran is quite dangerous and troublesome, because the existence of intrinsic “IEEE_xxx” modules etc apparently implies that NaN etc are treated “correctly” (at least for their test functions), but in fact, they are not (depending on compiler options, whose behavior cannot be predicted in a portable way).

1 Like

It will depend on the field you work in, but I would consider any code that “fails” when compiled with -Ofast or -ffast-math type options to be incorrect.

1 Like

For this code,

program main
    use ieee_arithmetic
    implicit none
    real :: x, y

    x = ieee_value( 1.0, ieee_quiet_nan )
    print *, x, sin( x )
    print *, isnan( x ), ieee_is_nan( x )

    y = ieee_value( 1.0, ieee_signaling_nan )
    print *, y, sin( y )
    print *, isnan( y ), ieee_is_nan( y )
end

ifort-2021.10 and ifx-2023.2 give with -standard-semantics + -O3 or -Ofast:

 NaN NaN
 T T
 NaN NaN
 T T

while gfortran (10,11,12,13) gives

(with -O3)
              NaN              NaN
 T T
              NaN              NaN
 T T
(with -Ofast or -ffast-math)
              NaN              NaN
 F F
              NaN              NaN
 F F

I think this is not something that depends on one’s work, but more related to the “optimization” that should not be made (IMHO).

To understand the behavior more, I’ve tried compiling this wrapper code:

!! mylib.f90
module mylib
    use ieee_arithmetic, only: ieee_is_nan
    implicit none
contains

function my_isnan( x ) result( res )
    real, intent(in) :: x
    logical :: res
    res = ieee_is_nan( x )
end
end module

Then, gfortran-12 -O3 gives this intermediate representation:

$ gfortran-12 -c -O3 -fdump-tree-original mylib.f90
$ cat mylib.f90.005t.original 

__attribute__((fn spec (". r ")))
logical(kind=4) my_isnan (real(kind=4) & restrict x)
{
  logical(kind=4) res;

  res = SAVE_EXPR <*x> unord SAVE_EXPR <*x>;
  return res;
}

which references the dummy argument x, while compilation with -ffast-math gives

$ gfortran-12 -c -ffast-math -fdump-tree-original mylib.f90
$ cat mylib.f90.005t.original

__attribute__((fn spec (". r ")))
logical(kind=4) my_isnan (real(kind=4) & restrict x)
{
  logical(kind=4) res;

  res = 0;
  return res;
}

which replaces the result of ieee_is_nan() unconditionally by 0.

Also, if I use the above mylib.f90 together with the following main.f90

!! main.f90
program main
    use ieee_arithmetic
    use mylib, only: my_isnan
    implicit none
    real :: x, y

    x = ieee_value( 1.0, ieee_quiet_nan )
    print *, x, sin( x )
    print *, isnan( x ), ieee_is_nan( x ), my_isnan( x )

    y = ieee_value( 1.0, ieee_signaling_nan )
    print *, y, sin( y )
    print *, isnan( y ), ieee_is_nan( y ), my_isnan( y )
end

simultaneous compilation of mylib.f90 and main.f90 gives a wrong result with -ffast-math:

$ gfortran-12 -O3 mylib.f90 main.f90 && ./a.out
              NaN              NaN
 T T T
              NaN              NaN
 T T T
$ gfortran-12 -ffast-math mylib.f90 main.f90 && ./a.out
              NaN              NaN
 F F F
              NaN              NaN
 F F F

while separate compilation with -O3 for mylib.f90 and -ffast-math for main.f90 gives T only for my_isnan():

$ gfortran-12 -c -O3 mylib.f90
$ gfortran-12 -ffast-math mylib.o main.f90 && ./a.out
              NaN              NaN
 F F T
              NaN              NaN
 F F T

The description says it all:

-ffast-math

Sets the options -fno-math-errno, -funsafe-math-optimizations, -ffinite-math-only, -fno-rounding-math, -fno-signaling-nans, -fcx-limited-range and -fexcess-precision=fast.

This option causes the preprocessor macro __FAST_MATH__ to be defined.

This option is not turned on by any -O option besides -Ofast since it can result in incorrect output for programs that depend on an exact implementation of IEEE or ISO rules/specifications for math functions. It may, however, yield faster code for programs that do not require the guarantees of these specifications.

In particular, the sub-flag:

-ffinite-math-only

Allow optimizations for floating-point arithmetic that assume that arguments and results are not NaNs or ±Infs.

So IMO, what you’d like to achieve is fundamentally incompatible with this compiler flag.

My own understanding has been that -Ofast is to be used only on problems which are numerically well-conditioned, suffer no harm from rearranging expressions (-Ofast also turns on -fno-protect-parens, meaning parentheses are not strict) and for which you assume responsibility yourself there will be no NaN or Inf.

In other words:

  • compiling the is_nan, is_finite, is_inf, is_posinf routines with -ffast-math, you may as well set the output to a constant .false., .true., .false., .false. (use can use the __FAST_MATH__ preprocessor flag)
  • compiling a numerical kernel with -ffast-math, and then performing a Inf/NaN test (not compiled with -Ofast) on the output results, also makes no sense; you just assured the compiler everything will be finite in the numerical kernel, so why would you want to test it

Unfortunately, the advice given on the web on -Ofast is similar to this sign located at The Storr on the Isle of Skye, Scotland:


(Source: Old Man of Storr - Walk Guide & Photography — And Then I Met Yoko | Visual Travel Guides)

Actually, what I feel strongly problematic is more about the description of this kind of things in the docs and man pages, for example,

https://www.intel.com/content/www/us/en/docs/fortran-compiler/developer-guide-reference/2023-0/ieee-is-nan.html

This man page says that

Elemental Module Intrinsic Function (Generic): Returns whether an IEEE value is Not-a-Number (NaN).

and never says that “it may not behave as you expect” caveat things. I assumed that if a compiler is standard-compliant / conforming, it follows the above kind of description. Is it not correct? (I.e., do I need to check every options to see whether some function is guaranteed to work in a portable way or not, if compilers are standard-conforming?)

In other words this option basically says “I am no longer a standards conforming compiler at this point.” (at least in some specific ways)

so

For Intel’s compiler documentation this may be a valid point, but they are basically just copying from the standard here. The standard already assumes standards conformance. I.e. the whole document has an implicit “things may not behave as expected if your compiler doesn’t conform to the standard.” So if you tell your compiler “it’s ok if you don’t strictly adhere to the standard” it’s kind of expected to get some unexpected behavior. Not to mention the fact that the code you’re trying to write/use (i.e. is_nan) already doesn’t serve much purpose if the rest of the program stays within the bounds of the standard.

1 Like

To be fair, -Ofast and -ffast-math both say that they do non-standard optimizations, like rearranging arithmetic, etc.

I understood @septc’s answer was targetting the Intel Classic compiler - ifort, where the default settings are “aggressive”, and -standard-semantics are NOT turned on by default.

The new ifx compiler has changed some of the default options, in particular the default floating point behaviour. An overview can be found in this document: Porting Guide for Intel® Fortran Compiler.

1 Like

As mentioned, the “why” is not what I am interested in. But since many ask, let me explain a bit.

  1. Why would I like to test Inf and NaN correctly even when compilers are invoked with aggressive optimization flags?

This is the answer:

I hold exactly the same opinion. I do not code for myself. I code for potential users whom I do not know. I cannot prevent them from compiling my code with aggressive optimization flags. It is my responsibility to ensure that my code works correctly regardless of the compiler flags. I do not think this is a very bizarre situation for those who believe that they are developing software rather than coding for themselves (there is a difference between the two).

  1. Does -Ofast do any good in any situation?

I do not know in general, but I do know it is not beneficial in the applications that my code targets. However, I cannot prevent my users from applying this flag. There must be a reason for the existence of such flags.

  1. Is it impossible to test Inf and NaN correctly under aggressive optimization flags (fundamentally incompatible with such compiler flags)?

I do not think so. As mentioned, I have some functions that worked well until gfortran 13. The same functions work correctly with ifort, ifx, nagfor, Absoft af95, Classic Flang, AMD AOCC Flang, Huawei Bisheng Flang, Arm Flang, nvfortran, Oracle sunf95, and g95 (basically all solvers available to me except for gfortran-13), each of which is invoked with the most aggressive optimization flag that I can conceive according to their documentations.

1 Like

You could inspect the result of compiler_options and issue a warning or error.

If you can’t trust the users to read the build documentstion, you could distribute it in binary form.

In Python, one can also do things that are wrong, but there is a common (and voluntary) rule: Code Style — The Hitchhiker's Guide to Python

Thanks for your advice.

I prefer to give users the most flexibility and assume myself the responsibility of ensuring the code works correctly.

Anyway, when developing an open-source package, I believe it is a good idea to remember that others may reuse (part of) your code without your building system.

Does the solution provided by @RonShepard in the other thread work? Is this expected? - #5 by RonShepard

Well, the other thread is on a quite different problem. The suggested solution is nice, but does not help with the problem under discussion.