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.
gfortran optimisation of overflowing integers
-
15-06-2023 - |
문제
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?
해결책
다른 팁
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.
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