Geek Answers Handbook

Infinity in Fortran

What is the safest way to set the + Infinity variable in Fortran? At the moment I am using:

program test implicit none print *,infinity() contains real function infinity() implicit none real :: x x = huge(1.) infinity = x + x end function infinity end program test

but I wonder if there is a better way?

+6
fortran infinity
source share
5 answers

If your compiler supports ISO TR 15580 IEEE Arithmetic, which is part of the so-called Fortran 2003 standard, than you can use the procedures from ieee_ * modules.

 PROGRAM main USE ieee_arithmetic IMPLICIT NONE REAL :: r IF (ieee_support_inf(r)) THEN r = ieee_value(r, ieee_negative_inf) END IF PRINT *, r END PROGRAM main
+7
source share

I'm not sure if this solution works on all compilers, but this is a good mathematical way to achieve infinity like -log (0).

 program test implicit none print *,infinity() contains real function infinity() implicit none real :: x x = 0 infinity=-log(x) end function infinity end program test

Also works well for complex variables.

+1
source share

I do not know about security, but I can offer you an alternative method. I learned this way:

 PROGRAM infinity IMPLICIT NONE INTEGER :: inf REAL :: infi EQUIVALENCE (inf,infi) !Stores two variable at the same address DATA inf/z'7f800000'/ !Hex for +Infinity WRITE(*,*)infi END PROGRAM infinity

If you use exceptional values ​​in expressions (I don’t think it is recommended at all), you should pay close attention to how your compiler deals with them, otherwise you may get unexpected results.

0
source share

I would not rely on the compiler to support the IEEE standard, and do what you did with two changes:

  • I would not add huge(1.)+huge(1.) , Since on some compilers you can end with -huge(1.)+1 ---, and this can lead to a memory leak (I don’t know the reason, but this is an experimental fact, so to speak).

  • Here you are using real types. I personally prefer to store all of my floating point numbers as real*8 , so all floating point constants have the value d0 , for example: huge(1.d0) . Of course, this is not the rule. some people prefer to use both real -s and real*8 -s.

0
source share

It seems to work for me. Define a parameter

 double precision,parameter :: inf = 1.d0/0.d0

Then use it in tests.

  real :: sng double precision :: dbl1,dbl2 sng = 1.0/0.0 dbl1 = 1.d0/0.d0 dbl2 = -log(0.d0) if(sng == inf) write(*,*)"sng = inf" if(dbl1 == inf) write(*,*)"dbl1 = inf" if(dbl2 == inf) write(*,*)"dbl2 = inf" read(*,*)

When compiling with ifort and running I get

 sng = inf dbl1 = inf dbl2 = inf
0
source share

More articles:


All Articles