gfortran optimisation of overflowing integers

Question

I have "inherited" a simple congruential pseudo-random-number generator that looks like this:

  subroutine ribm(rndm,ial)
    implicit none
    integer :: ial
    real :: rndm, al
    ial=ial*65539
    if(ial.lt.0) then
      ial=ial+2147483647+1
    endif
    al=ial
    rndm=al*0.4656613e-9
  end subroutine ribm

The ial variable contains the RNG seed and it is initialised once outside the generator. The generator is supposed to produce pseudo-random numbers between 0 and 1, as usual. Therefore, ial must always stay non-negative. The code tries to enforce this condition using the if statement.

I test the generator using the following program:

program testribm
  implicit none
  real :: r
  integer :: i
  data i /12345/

  call ribm(r,i)
  call ribm(r,i)
  call ribm(r,i)
  call ribm(r,i)
  call ribm(r,i)
  call ribm(r,i)
  print *,r,i
contains
  subroutine ribm(rndm,ial)
  ! [...]
  end subroutine ribm
end program testribm

If I compile with gfortran v4.8.2 with optimisation -O1 or lower, I get:

 0.580895185      1247462939

OTOH, if I use -O2 I get:

-0.709615409     -1523887535

Note the negative seed value. Optimisation affects the integer arithmetics and the final result of the repeated calls to ribm.

I know how to work around/fix this problem. My question is: is this a gfortran bug? Or is this behaviour due to the fact that the code is non-standard? But if this is the case, shouldn't gfortran produce some warning message?

Answer 1

Optimization can change the actual calculations done. Don't be surprised if the results change for a non-standard calculation, such as one that involves integer overflow. My suggestion to implement algorithms in Fortran that use unsigned integers is to use the next-larger integer type. e.g., if the algorithm should be done with unsigned 32-bit integers, use Fortran's signed 64-bit integer.

Answer 2

I would add to M.S.B.'s answer, that the problem occurs only when the compiler optimizes several calls to ribm in one go. (My previous suggestion was invalid.)

But the answer to the original question. I don't think it's a compiler bug. 2147483647+1 is simply out of range. The question is why there is no warning even with -pedantic -Wall -Wextra. It turns out, pointed out by @Arek' Fu that one needs -Wstrict-overflow, whis is included only as -Wstrict-overflow=1 in -Wall.

Where possible, I recommend using functions iand, ior, ieor, ishft and similar for portable results, instead of arithmetic.

Answer 3

both answers are nice, as for a solution, I would think using parenthesis is the right way to go to prevent the compiler optimizing away?

subroutine ribm(rndm,ial)
implicit none
integer :: ial
real :: rndm, al
ial=ial*65539
if(ial.lt.0) then
  ial=(ial+1)+2147483647
endif
al=ial
rndm=al*0.4656613e-9
end subroutine ribm