Infinity in Fortran

Question

What is the safest way to set a variable to +Infinity 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 am wondering if there is a better way?

Solution 1

If your compiler supports ISO TR 15580 IEEE Arithmetic which is a part of so-called Fortran 2003 standard than you can use 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

Solution 2

I'm not sure if the solution bellow works on all compilers, but it's a nice mathematical way of reaching infinity as -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 nicely for complex variables.

Solution 3

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

1. I would not add huge(1.)+huge(1.), since on some compilers you may end up with -huge(1.)+1 --- and this may cause a memory leak (don't know the reason, but it is an experimental fact, so to say).

2. You are using real types here. I personally prefer to keep all my floating-point numbers as real*8, hence all float constants are qualified with d0, like this: huge(1.d0). This is not a rule, of course; some people prefer using both real-s and real*8-s.

Comments

[Vladimir F Героям слава]

real*8 and double precision (1.d0) are not necesarilly the same real kind. And of course, whether to use single or double precision is not a question of personal preference, but of mathematical arguments and tests.

[Vladimir F Героям слава]

It will work on IEEE compliant machines, if you do not enable FPE trapping. But if you work with infinities, it would be silly to do so.

[Paul Wintz]

This does not work, at least in gfortran. It raises a compiler error.