Infinity in Fortran

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开发者_运维知识库?

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

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.

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.

I don't know about safest, but I can offer you an alternative method. I learned to do it 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 are using exceptional values in expressions (I don't think this is generally advisable) you should pay careful attention to how your compiler handles them, you might get some unexpected results otherwise.

This seems to work for me. Define a parameter

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

Then use it in if 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 compiled with ifort & run, I get

sng = inf
dbl1 = inf
dbl2 = inf