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?

Subcontinent

Answer 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

Revivalist

Answer 2

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.

Micropathology

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. – Mercenary

Answer 3

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.

Novgorod

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. – Mercenary

Answer 4

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.

Ministration

Answer 5

For anyone coming to this thread, the hexadecimal conversion to define inf/-inf can be easily done by:

use iso_fortran_env, only: sp => real32, dp => real64, qp => real128
!> single precision infinities
real(sp) :: pinf_sp = real(z'7f800000',sp) 
real(sp) :: ninf_sp = real(z'ff800000',sp) 
!> double precision infinities
real(dp) :: pinf_dp = real(z'7ff0000000000000',dp)
real(dp) :: ninf_dp = real(z'fff0000000000000',dp)
!> quadruple precision infinities
real(qp) :: pinf_qp = real(z'7fff0000000000000000000000000000',qp)
real(qp) :: ninf_qp = real(z'ffff0000000000000000000000000000',qp)

European

Answer 6

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

Ott

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