Is there a standard way to check for Infinite and NaN in Fortran 90/95?

Question

I've been trying to find a standards-compliant way to check for Infinite and NaN values in Fortran 90/95 but it proved harder than I thought.

• I tried to manually create Inf and NaN variables using the binary representation described in IEEE 754, but I found no such functionality.
• I am aware of the intrinsic ieee_arithmetic module in Fortran 2003 with the ieee_is_nan() and ieee_is_finite() intrinsic functions. However it's not supported by all the compilers (notably gfortran as of version 4.9).

Defining infinity and NaN at the beginning like pinf = 1. / 0 and nan = 0. / 0 seems hackish to me and IMHO can raise some building problems - for example if some compilers check this in compile time one would have to provide a special flag.

Is there a way I can implement in standard Fortran 90/95?

function isinf(x)
! Returns .true. if x is infinity, .false. otherwise
...
end function isinf

and isnan()?

Lactoprotein

gnu fortran 4.10 fixes this – Debenture

GCC 5 and more recent do support IEEE_ARITHMETIC, but support of older versions is still an issue and will continue to be for a long time. – Encrimson

Answer 1

The simple way without using the ieee_arithmatic is to do the following.

Infinity: Define your variable infinity = HUGE(dbl_prec_var) (or, if you have it, a quad precision variable). Then you can simply check to see if your variable is infinity by if(my_var > infinity).

NAN: This is even easier. By definition, NAN is not equal to anything, even itself. Simply compare the variable to itself: if(my_var /= my_var).

Indicate

Odd that nobody else mentioned the NaN case. I guess they were worried that you might be using fortran on a non-ieee cpu? Aside from that, the NaN check is very robust. The infinity check, on the other hand, is a bit less elegant, since it depends on the data type. – Podite

@amaurea: You can take care of the data-type issue by using INTERFACE and linking the few types you might be using into the same MODULE PROCEDURE. – Indicate

@IanH: You could reasonably set infinity = 1.e100_dp and still obtain a good infinity checker. Even when using small scales (like 1e-24 m) to describe big scales (like 1e+24 m), you still are no where near 1e100. – Indicate

Thank you for the answer, and good catch about the MODULE PROCEDURE. One can make a very clean interface with it. If my code is worth it I might release it. Thanks again! – Lactoprotein

An optimizer will indeed optimize the check out if fast_math or similar is requested. Alas, it will even optimize out the gfortran's isnan intrinsic. It is the programmer's responsibility to take care about this when he requests unsecure fast math and still wants to detect NaN. – Encrimson

HUGE is the largest non-infinity; the condition myvar > HUGE(myvar) should only be true if myvar is infinite, but naming HUGE(myvar) as infinity is misleading. – Margrettmarguerie

How much of this answer is functional vs. standard-compliant? In other words, is this how it ought to be done, or is this just how it can be done? – Perl

@Perl using iee_arithmetic is how it should be done (so if your compiler doesn't support it, get a different/better compiler). This is a standards-compliant way that is reasonably functional: variables are "infinite" when larger than a value that doesn't make sense for the code (and infinitesimal when smaller than a certain value); the NaN check could be optimized with some compilers though. – Indicate

I tried to check for NaN in nvfortran using if(my_var /= my_var), but it didn't work. Hopefully, the alternative expression if((my_var-my_var) /= 0) did work in the newer compilers (2022/2023). NOTE: I know all is supposed to work, but I can provide a code example showing that if(my_var /= my_var) fails. – Mysia

Answer 2

I don't have enough rep to comment so I'll "answer" regarding Rick Thompson's suggestion for testing infinity.

if (A-1 .eq. A) 

This will also be true if A is a very large floating point number, and 1 is below the precision of A.

A simple test:

subroutine test_inf_1(A)
    real, intent(in) :: A
    print*, "Test (A-1 == A)"
    if (A-1 .eq. A) then
        print*, "    INFINITY!!!"
    else
        print*, "    NOT infinite"
    endif
end subroutine

subroutine test_inf_2(A)
    real, intent(in) :: A
    print*, "Test (A > HUGE(A))"
    if (A > HUGE(A)) then
        print*, "    INFINITY!!!"
    else
        print*, "    NOT infinite"
    endif
end subroutine


program test
    real :: A,B

    A=10
    print*, "A = ",A
    call test_inf_1(A)
    call test_inf_2(A)
    print*, ""

    A=1e20
    print*, "A = ",A
    call test_inf_1(A)
    call test_inf_2(A)
    print*, ""

    B=0.0 ! B is necessary to trick gfortran into compiling this
    A=1/B
    print*, "A = ",A
    call test_inf_1(A)
    call test_inf_2(A)
    print*, ""

end program test

outputs:

A =    10.0000000    
Test (A-1 == A)
    NOT infinite
Test (A > HUGE(A))
    NOT infinite

A =    1.00000002E+20
Test (A-1 == A)
    INFINITY!!!
Test (A > HUGE(A))
    NOT infinite

A =          Infinity
Test (A-1 == A)
    INFINITY!!!
Test (A > HUGE(A))
    INFINITY!!!

Austenite

Answer 3

I have used:

  PROGRAM MYTEST
  USE, INTRINSIC :: IEEE_ARITHMETIC, ONLY: IEEE_IS_FINITE      
  DOUBLE PRECISION :: number, test
  number = 'the expression to test'
  test = number/number
  IF (IEEE_IS_FINITE(test)) THEN
     WRITE(*,*) 'We are OK'
  ELSE
     WRITE(*,*) 'Got a problem'
  END IF         
     WRITE(*,*) number, test
  END PROGRAM MYTEST

This will print 'Got a problem' for number = 0.0D0, 1.0D0/0.0D0, 0.0D0/0.0D0, SQRT(-2.0D0), and also for overflows and underflows such as number = EXP(1.0D800) or number = EXP(-1.0D800). Notice that generally, things like number = EXP(1.0D-800) will just set number = 1.0 and produce a warning at compilation time, but the program will print 'We are OK', which I find acceptable.

OL.

Mallett

The question asks for NotANumber, but your exampletests for finiteness. It will also catch +Inf and -Inf as false positives. Also notice that the question was really mainly interested about what to do when IEEE_ARITHMETIC is not available and asks how to do that using Fortran 90/95. – Encrimson

In addition to Vladimir F's concerns, there's no guarantee that even with ieee_arithmetic that ieee_support_datatype(test) is true. If it isn't, it is not permitted to consider ieee_is_finite(test). – Vinasse

I was bit too harsh, the question does also ask for an equivalent of IEEE_IS_FINITE(). But the point is that the OP knows it exists but asks for an alternative. – Encrimson

Answer 4

No.

The salient parts of IEEE_ARITHMETIC for generating/checking for NaN's are easy enough to write for gfortran for a particular architecture.

Liaoyang

About your second statement, they don't seem to agree gcc.gnu.org/ml/gcc-bugs/2012-10/msg00580.html or am I missing something? – Lactoprotein

Btw I added a link to the relevant gfortran bug (which dates from 2006 and has status NEW). – Lactoprotein

Writing all of IEEE_ARITHMETIC for all architectures that gfortran supports, which is what that bug deals with, would be difficult! Writing the select bits that generate/check for NaN's, on a particular architecture (for example, x64) is pretty easy. See sites.google.com/site/tprincesite/Home/gfortran-ieee-arithmetic for one example by Tim Prince that relies on the inequality of any two NaN's, an alternative approach is to use the Fortran bit intrinsics to generate/test for the specific patterns that indicate a NaN. (See the "original" gfortran bug 29383 for more discussion.) – Liaoyang

Answer 5

for testing NaN none of the things worked eg.if testing real s2p to see if it is NaN then

if(isnan(s2p)) 

did not work in gfortran nor did

if(s2p.ne.s2p). 

The only thing that worked was

if(.not.s2p<1.and..not.s2p>1)  

though to make real sure u may want to add

if(.not.s2p<1.and..not.s2p>1.and..not.s2p==1)    

Lentil

Answer 6

No.

Neither is there a standards-compliant way of checking for infinities or NaNs in Fortran 90/95, nor can there be a standards-compliant way. There is no standards-compliant way of defining either of these quasi-numbers in Fortran 90/95.

Isiah

Well, in the scope of the IEEE Floating Point standard, that is. – Lactoprotein

Answer 7

For Fortran, 1/infinity=0 thus, divide your variable by zero i.e

program test
implicit none
real :: res
integer :: i

do i=1,1000000
    res=-log((1.+(10**(-real(i))))-1.)
    print*,i,res
    if ((1./res)==0.) then
        exit
    end if
end do

end program

there's your infinity check. No complication necessary.

Pedicab

Answer 8

For Inf it seems to work that if (A-1 .eq. A) is true, then A is Inf

Redfaced

Maybe you might want to elaborate a bit more. – Margerymarget