Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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

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()?

like image 770
astrojuanlu Avatar asked Jun 30 '13 11:06

astrojuanlu


People also ask

How do you define NaN in Fortran?

Not a Number (NaN) is represented by the largest value that the exponent can assume (all ones), and a nonzero fraction. Normalized REAL and DOUBLE PRECISION numbers have an implicit leading bit that provides one more bit of precision than is stored in memory.


2 Answers

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

like image 94
Kyle Kanos Avatar answered Sep 18 '22 18:09

Kyle Kanos


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!!! 
like image 23
Ethan Gutmann Avatar answered Sep 21 '22 18:09

Ethan Gutmann