I am puzzled again by the behavior of ifort, which (sometimes) evaluates 1 / Inf to NaN. I thought it should always be zero.
Given the experience with another question, I think it might be better to ask the experts here before calling it a bug.
Code:
! test.f90
program test
use, intrinsic :: ieee_arithmetic, only : ieee_value, ieee_positive_inf, ieee_support_inf
implicit none
real :: a(5)
real :: b
a = 1.0
b = ieee_value(b, ieee_positive_inf)
write (*, *) 'IEEE_SUPPORT_INF = ', ieee_support_inf(a)
write (*, *) 'A = ', a
write (*, *) 'B = ', b
write (*, *) 'A / B = ', a / b
end program test
Result:
$ uname -a && ifort --version && ifort test.f90 && ./a.out
Linux 5.15.0-52-generic #58-Ubuntu SMP Thu Oct 13 08:03:55 UTC 2022 x86_64 x86_64 x86_64 GNU/Linux
ifort (IFORT) 2021.7.1 20221019
Copyright (C) 1985-2022 Intel Corporation. All rights reserved.
IEEE_SUPPORT_INF = T
A = 1.000000 1.000000 1.000000 1.000000 1.000000
B = Infinity
A / B = NaN NaN NaN NaN 0.0000000E+00
See my experiment with Compiler Explorer. Strangely, I cannot reproduce the results above without specifying the -O option. On my computer, no option was needed.
What do you think is going on here? According to previous discussions, thanks to @greenrongreen , we understand why the fifth entry of A / B differs from others, but why are the first four NaN?
Do Fortran standards specify what 1 / Inf should be, or are compilers free to (un)define it?
From your point of view, is this (default) behavior acceptable / reasonable? What would be your choice if you were going to make a compiler?
You need to use -fp-model=strict to get the semantics you want.
D:\Projects>ifort /fp:strict t.f90
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.7.1 Build 20221019_000000
Copyright (C) 1985-2022 Intel Corporation. All rights reserved.
Microsoft (R) Incremental Linker Version 14.33.31630.0
Copyright (C) Microsoft Corporation. All rights reserved.
-out:t.exe
-subsystem:console
t.obj
D:\Projects>t.exe
IEEE_SUPPORT_INF = T
A = 1.000000 1.000000 1.000000 1.000000 1.000000
B = Infinity
A / B = 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00
0.0000000E+00
While most compilers would set 1 / inf = 0 by default, do we consider producing 1 / inf = NaN an acceptable & reasonable (default) behavior rather than something to worry about? If the answer is a firm yes, I will just take that note and learn to live with it, even though this is the first time I ever hear or even think about it.
Optimizing compilers orften have default behaviors that favor performance over edge-cases, such as full IEEE semantics. The Fortran standard is completely silent on this topic. In the case of ifort, you need the strict option if you will be using anything from IEEE_EXCEPTIONS and possibly even IEEE_ARITHMETIC.
How realistic is it to expect a real-world application to be doing 1/Inf?
If the infinity results from a FP overflow, like in A/(B*C), I would say it can be expected and it would make sense to get 0 from this operation and continue. Consider for example:
real :: a=1e-15, b=1e35, c=1e5
print *, a/b ! obviously 0.0
print *, a/(b*c) ! I'd expect 0.0 again, not a NaN
Edit:
Worse enough, according to the Standard, a/b/c can be computed by any mathematically equivalent expression (in the absence of parentheses), so either as (a/b)/c - surely 0.0 or a/(b*c) - now as we see 0.0 or NaN. Not nice.
For that value of a, maybe 0.0 is the “best” result. But for large values of a, say anywhere in the range [1e5,1e35], 0.0 would be an unexpected value to return. The compiler is not required to assess every possible way to evaluate the result in order to avoid overflows or underflows, it is allowed to just pick a way (that respects any parentheses constraints) and evaluate it that way. Returning NaN at least tells the programmer that its chosen evaluation method resulted in something unexpected. The programmer, with knowledge of the meaning and ranges of the variables, can then add parentheses or other alternative expressions as appropriate.
@zaikunzhang , it is unclear if there is much anyone can “teach” you here. The issues you are dealing with are in the processor-dependent category and chances are very high there will considerable divergence in views on the exact approaches to follow even if there is some consensus on the broad “brush strokes”. This is especially going to be the case among compiler implementations and their developers with floating-point operations. You will not get Intel and/or gfortran or another compiler to all do the same things by default.
Besides your immediate concerns are specific to Intel Fortran. You will notice the good folks at Intel are effectively informing you to decide on the best course of action for yourself - whether to use aggressive optimizations per Intel’s implementation with its default of -fp-model=fast option or go with some other options that can better suit your needs. They would want you to study their documentation in-depth, also try things out, and decide for yourself what options suit your needs: the discourse is limited.