gfortran: why "-posinf" leads to Arithmetic overflow?

Question

to get +-inf on 64 bit system i used the next code

double precision, parameter :: pinf           = transfer(z'7FF0000000000000',1d0) ! 64 bit
double precision, parameter :: ninf           = transfer(z'FFF0000000000000',1d0) ! 64 bit

and it works well.

On 32-bit I've got an compilation error only(!) for ninf:

double precision, parameter :: ninf           = transfer(z'FFF0000000000000',1d0
                                                                           1     
Error: Integer too big for integer kind 8 at (1)

assignment ninf = -pinf not helps and leads to compilation Arithmetic overflow error:

double precision, parameter :: ninf           = -pinf
                                                         1
Error: Arithmetic overflow at (1)

I know about ieee_arithmetic module but gcc don't handle it.

Is there any multi-architecture way to set constants to positive/negative infinities?

Update Gfortran option -fno-range-check suppress errors and successfully compile that code.

It's not important but I'm still interesting. Why gfortran allows constant definition of +Infinity but yelling in loud about exactly the same thing with -Infinity?

Answer 1

In this case gfortran is internally representing your hexadecimal ("Z") literals as the largest unsigned integer size available. Since transfer is a Fortran intrinsic, and Fortran does not have unsigned integers, the first thing gfortran does is to assign the literal to a signed type, which causes your bit pattern for negative infinity to overflow. This happens in many other cases where you use BOZ literals, and I think that this is a bug in gfortran.

I think this only shows up on a 32 bit system because on your 64 bit system, gfortran probably has a 128 bit integer type available; a 128 bit signed integer will not "overflow" with that bit pattern.

But it is also the case that your code does not conform to the Fortran standard, which says that hex literals can only appear inside data statements or the functions int, real, or dble. However, putting a hex literal in dble does the same thing as transfer anyway. If gfortran did not have a bug in it, your program would work, but it would technically be incorrect.

Anyway, the following code works for me in gfortran, and I believe it will solve your problem in a way that's standard-compliant and avoids -fno-range-check:

integer, parameter :: i8 = selected_int_kind(13)
integer, parameter :: r8 = selected_real_kind(12)
integer(i8), parameter :: foo = int(Z'7FF0000000000000',i8)
integer(i8), parameter :: bar = ibset(foo,bit_size(foo)-1)
real(r8), parameter :: posinf = transfer(foo,1._r8)
real(r8), parameter :: neginf = transfer(bar,1._r8)
print *, foo, bar
print *, posinf, neginf
end

Output:

  9218868437227405312    -4503599627370496
                  Infinity                 -Infinity

The key is to create the pattern for positive infinity first (since it works), and then create the pattern for negative infinity by simply setting the sign bit (the last one). The ibset intrinsic is only for integers, so you then have to use transfer on those integers to set your real positive/negative infinity.

(My use of i8/r8 is just habit, since I've worked with compilers where the kind parameter was not equal to the number of bytes. They are both equal to 8 in this case.)

Answer 2

I'm not using the same compiler as you are (I'm using g95 with compiler option -i4 set for 32-bit integers, and one workaround (if you're staunch about using transfer for that purpose) that I found was to specify the integer argument as a parameter like so:

Note: with my compiler, I was able to assign the number directly to the parameter. I'm not sure if it's the same on yours, but I'm pretty sure that you're only really supposed to use the transfer function when you're not really dealing with constants -- like if you're doing fancy stuff with floating point numbers and need like really nitty gritty control over the representation thereof.

Note the variables pdirect and ndirect.

program main
integer(8), parameter :: pinfx= z'7FF0000000000000'
integer(8), parameter :: ninfx= z'FFF0000000000000'
double precision, parameter :: pinf = transfer(pinfx, 1d0)
double precision, parameter :: ninf = transfer(ninfx, 1d0)

double precision, parameter :: pdirect = z'7FF0000000000000'
double precision, parameter :: ndirect = z'7FF0000000000000'


write (*,*) 'PINFX  ', pinfx
write (*,*) 'NINFX  ', ninfx
write (*,*) 'PINF   ', pinf
write (*,*) 'NINF   ', ninf
write (*,*) 'PDIRECT', pdirect
write (*,*) 'NDIRECT', ndirect

end program

This produces the output:

 PINFX   9218868437227405312
 NINFX   -4503599627370496
 PINF    +Inf
 NINF    -Inf
 PDIRECT +Inf
 NDIRECT +Inf

I hope this helps!