Dear all,
This may be really a lazy question, I need some function or subroutine to tell me if the variable, say X, is NaN, +infinity or -infinity.
I know there is an intrinsic function isnan
isnan(X)
which return True of False therefore tell me if X is NaN. .
But how about +infinity or -infinity? Is there some similar intrinsic function or trick can be used?
something like
isinfinity(X) ?
Thank you very much in advance!
You should use the IEEE_ARITHMETIC intrinsic module. It defines all manner of functions and data types for just this sort of things.
isnan() is a non-standard function. The standard one is ieee_is_nan() . To check for finite/infinite numbers you have ieee_is_finite(). And there are many others.
The intrinsic module ieee_arithmetic contains procedures to do these jobs.
However, note that they may not always work as expected when the compilers are invoked with some “special” options. For example, see the following discussions.
This program defines functions is_plus_infinity and is_minus_infinity
module infinity_mod
use, intrinsic :: ieee_arithmetic
implicit none
integer, parameter :: wp = kind(1.0d0)
contains
elemental function is_plus_infinity(x) result(tf)
real(kind=wp), intent(in) :: x
logical :: tf
tf = (.not. ieee_is_finite(x)) .and. (.not. ieee_is_negative(x))
end function is_plus_infinity
!
elemental function is_minus_infinity(x) result(tf)
real(kind=wp), intent(in) :: x
logical :: tf
tf = (.not. ieee_is_finite(x)) .and. ieee_is_negative(x)
end function is_minus_infinity
end module infinity_mod
program test_infinity
use infinity_mod
implicit none
real(kind=wp) :: z, vec(2)
z = 0.0_wp
vec = [-1.0_wp,1.0_wp]/z
print "(a15,*(f15.12))","value",vec
print "(a15,*(l15))","+Inf?",is_plus_infinity(vec)
print "(a15,*(l15))","-Inf?",is_minus_infinity(vec)
end program test_infinity
and gives output
value -Infinity Infinity
+Inf? F T
-Inf? T F
Compiling the code with gfortran -Ofast with gfortran 9.3.0 on Ubuntu 20.04, I got
value -Infinity Infinity
+Inf? F F
-Inf? F F
This is not surprising due to the aggressive optimization of -Ofast.
NaN are fun: a NaN value is the only value not equal to itself! As it is not a value… Therefore you can use a simple test:
real :: x
x = -1
x = sqrt(x)
print *, x ! Should print NaN
if (x /= x) print *, "Not a Number!"
In terms of coding practice, my experience is that testing for NaN and infinity is usually the wrong thing to do. In the past, I would do some “risky” computation that might produce NaN or inf, and then test afterwards if I got one of these results and take some follow-up action. This almost always created more problems down the road, and the better approach was just to make the “risky” computation more robust. Sometimes that means improving the algorithm. Sometimes that means validating the input and stopping if it’s bad.
There are definitely valid reasons to use or test for NaN and inf, but my own experience is that these are the exception.
Thank you for the comments.
It was gfortran 9.3.0 on Ubuntu 20.04 (I tested just now also 10.3.0 on Ubuntu 20.04).
To be more precise,
$ gfortran --version && uname -a && lscpu | grep Nom
GNU Fortran (Ubuntu 9.3.0-17ubuntu1~20.04) 9.3.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions. There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
Linux 5.11.15-051115-generic #202104161034 SMP Fri Apr 16 10:40:30 UTC 2021 x86_64 x86_64 x86_64 GNU/Linux
Nom de modèle : Intel(R) Core(TM) i7-10610U CPU @ 1.80GHz
Thank you all! @Arjen @zaikunzhang @Beliavsky @vmagnin @nshaffer !
That is great!
If using IEEE_ARITHMETIC intrinsic module may lead to discrepancy due to -Ofast flag.
Then perhaps we can just use below to check NaN, +inf, -inf,
About if x is NaN,
like @vmagnin suggested, we can define a function to check if x = x or not.
About if x is +inf or -inf, can we use intrinsic HUGE function like below?
if x > huge(x), it means x is +inf
if x< -huge(x), it means is -inf
By the way, I agree with @nshaffer , when such checks are needed, the best solid fix is to check the code and use math to prevent such things from happening.
This will not work with gfortran -Ofast either — to be specific, gfortran 9.3.0 and 10.3.0 on Ubuntu 20.04 with Linux kernel 5.4.0-96-generic, the cpu being Intel(R) Core(TM) i7-4790 CPU @ 3.60GHz. However, not unlikely, it may be just because I have a flawed machine, bugged system, or simply bad luck. You may have a try on your side. Thanks.
Let me stress again that there is nothing wrong in this case — neither the code, nor gfortran. The outcome is expected if we read (and understand) the documentation as reminded and as we always should.
If you want some replacement that works (for the time being), you may try my GitHub repo mentioned above and in the other thread quoted before. They have been tested with 9 compilers (all that I could find) with aggressive optimization flags. I will be very happy to see where they do not work, and if possible, some suggested alternative implementation will be highly appreciated. Many thanks!
You should be able to simply write:
elemental function is_plus_infinity(x) result(tf)
real(kind=wp), intent(in) :: x
logical :: tf
tf = ieee_class(x) == IEEE_POSITIVE_INF
end function is_plus_infinity
elemental function is_minus_infinity(x) result(tf)
real(kind=wp), intent(in) :: x
logical :: tf
tf = ieee_class(x) == IEEE_NEGATIVE_INF
end function is_minus_infinity
Not quite! Just about every self-respecting optimiser will see that this can never be true, so it eliminates the condition.
Actually, I recently found out that the optimiser in gfortran (-O2) can eliminate an entire nested do-loop, rendering the program very fast indeed. The do-loop was intended to measure performance and I thought I could outsmart the optimiser by putting in some dummy calculation. It was correct in doing so and I was incorrect assuming I could fool it in that way.
I would not say “self-respecting” as the assumption that it can never be true is technically false… 
Could someone give an example where ieee_is_finite(x) and ieee_is_normal(x) give different results? Below is a code I wrote that uses some functions from the ieee_arithmetic intrinsic module.
program test_ieee_arithmetic
use, intrinsic :: ieee_arithmetic
implicit none
character (len=20) :: fmt = "(a16,*(l9))"
integer, parameter :: wp = kind(1.0) ! same output for wp = kind(1.0d0)
real(kind=wp) :: z,vec(6),NaN,t
z = 0.0_wp
t = tiny(z)
vec = [z/z,z/1.0_wp,-z/1.0_wp,1.0_wp/z,-1.0_wp/z,t**2]
print*,ieee_value(0.0,ieee_positive_inf),1.0/z ! 2 ways of geting +Inf
NaN = ieee_value(0.0,ieee_quiet_nan) ! get NaN
print*,"NaN == NaN?",NaN == NaN ! demonstrate that NaN /= NaN
print*
print "(19x,*(a9))","0.0/0.0","0.0/1.0","-0.0/1.0","1.0/0.0","-1.0/0.0","tiny^2"
print fmt,"ieee_is_nan" ,ieee_is_nan(vec)
print fmt,"ieee_is_negative",ieee_is_negative(vec)
print fmt,"ieee_is_finite" ,ieee_is_finite(vec)
print fmt,"ieee_is_normal" ,ieee_is_normal(vec)
end program test_ieee_arithmetic
! output:
! Infinity Infinity
! NaN == NaN? F
!
! 0.0/0.0 0.0/1.0 -0.0/1.0 1.0/0.0 -1.0/0.0 tiny^2
! ieee_is_nan T F F F F F
! ieee_is_negative F F T F T F
! ieee_is_finite F T T F F T
! ieee_is_normal F T T F F T
Here you go:
! denormal.f90 --
! Show a finite non-normal number
!
program denormal
use ieee_arithmetic
real :: x
x = tiny(x) / 10.0
write(*,*) x
write(*,*) 'Finite?', ieee_is_finite(x)
write(*,*) 'Normal?', ieee_is_normal(x)
end program denormal
Thanks. Generalizing the program to
program denormal
use ieee_arithmetic
implicit none
integer :: pow10
real :: x
write (*,"(a6,a10,3x,*(a12))") "pow10","x","is_finite","is_normal", &
"x==0.0","x<tiny","chk_normal"
do pow10 = 0,9
x = tiny(x) / (10.0**pow10)
write(*,"(4x,i1,2x,es10.4,*(l12))") pow10,x,ieee_is_finite(x), &
ieee_is_normal(x),x==0.0,x<tiny(x), (x>0.0 .and. x>=tiny(x)) .or. x == 0.0
end do
end program denormal
the output with gfortran is
pow10 x is_finite is_normal x==0.0 x<tiny chk_normal
0 1.1755E-38 T T F F T
1 1.1755E-39 T F F T F
2 1.1755E-40 T F F T F
3 1.1755E-41 T F F T F
4 1.1757E-42 T F F T F
5 1.1771E-43 T F F T F
6 1.1210E-44 T F F T F
7 1.4013E-45 T F F T F
8 0.0000E+00 T T T T T
9 0.0000E+00 T T T T T
and with Intel Fortran is
pow10 x is_finite is_normal x==0.0 x<tiny chk_normal
0 1.1755E-38 T T F F T
1 0.0000E+00 T T T T T
2 0.0000E+00 T T T T T
3 0.0000E+00 T T T T T
4 0.0000E+00 T T T T T
5 0.0000E+00 T T T T T
6 0.0000E+00 T T T T T
7 0.0000E+00 T T T T T
8 0.0000E+00 T T T T T
9 0.0000E+00 T T T T T
It appears that ieee_is_normal(x) gives the same result as
(x>0.0 .and. x>=tiny(x)) .or. x == 0.0
for nonnegative x. For Intel Fortran when x is less than tiny(x) it is set to zero, but for gfortran there are numbers where x is less than tiny(x) but not equal to zero.
Hm, what options do you use for Intel Fortran? Numbers x smaller than tiny(x) should be denormal, not zero. There is an option where such numbers are not gradually decreasing towards zero, but I would not expect that to be the default.
The output is
compiler version:
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.5.0 Build 20211109_000000
compiler options:
/nologo /F512000000 /o:ifort_denormal.exe
pow10 x is_finite is_normal x==0.0 x<tiny chk_normal
0 1.1755E-38 T T F F T
1 0.0000E+00 T T T T T
2 0.0000E+00 T T T T T
3 0.0000E+00 T T T T T
4 0.0000E+00 T T T T T
5 0.0000E+00 T T T T T
6 0.0000E+00 T T T T T
7 0.0000E+00 T T T T T
8 0.0000E+00 T T T T T
9 0.0000E+00 T T T T T
for
program denormal
use iso_fortran_env, only: compiler_version,compiler_options
use ieee_arithmetic
implicit none
integer :: pow10
real :: x
write (*,"('compiler version:',/,a)") trim(compiler_version())
write (*,"('compiler options:',/,a,/)") trim(compiler_options())
write (*,"(a6,a10,3x,*(a12))") "pow10","x","is_finite","is_normal", &
"x==0.0","x<tiny","chk_normal"
do pow10 = 0,9
x = tiny(x) / (10.0**pow10)
write(*,"(4x,i1,2x,es10.4,*(l12))") pow10,x,ieee_is_finite(x), &
ieee_is_normal(x),x==0.0,x<tiny(x), (x>0.0 .and. x>=tiny(x)) .or. x == 0.0
end do
end program denormal
I compiled your program with Intel oneAPI Fortran and the flag /Qftz- and got:
pow10 x is_finite is_normal x==0.0 x<tiny chk_normal
0 1.1755E-38 T T F F T
1 1.1755E-39 T F F T F
2 1.1755E-40 T F F T F
3 1.1755E-41 T F F T F
4 1.1757E-42 T F F T F
5 1.1771E-43 T F F T F
6 1.1210E-44 T F F T F
7 1.4013E-45 T F F T F
8 0.0000E+00 T T T T T
9 0.0000E+00 T T T T T