Bug with 32 bit compiler

I think that the following code reveals a bug with the 32 bit compiler. The variable a2 must always be positive.

module random_numbers
implicit none
integer, parameter :: dp = kind(1.d0)
real(dp), parameter :: pi = 4.d0*atan(1.d0)
contains

real(dp) function random_std_uniform()
! Function to generate a random number from the standard uniform distribution.
 do
   call random_number(random_std_uniform)
   if (random_std_uniform .gt. tiny(1.0)) exit
 end do
end function random_std_uniform

real(dp) function random_std_normal()
! Function to generate a random number from the standard normal distributon.  
  random_std_normal = sqrt(-2.d0*log(random_std_uniform())) * cos(2.d0*pi*random_std_uniform())
end function random_std_normal

end module random_numbers

program bug
use random_numbers
use iso_fortran_env
implicit none
integer i
real(dp) a, a2
write(*,*) 'Vern : ',compiler_version ()
write(*,*) 'Options : ',compiler_options ()
write(*,*)
! This do loop runs as expected.  a can be positive or negative, a**2 is always positive
write(*,*) 'OK with intermediate value'
do i = 1, 10
  a = random_std_normal()
  a2 = (a**2)
  write(*,*) i, a, a2
end do
write(*,*)
! This do loop which does not have an intermediate set does not run as expected with the 
! 32 bit compiler, it will on occasion return a negative number.
write(*,*) 'Bug appears without the intermediate value'
do i = 1, 10
  a2 = (random_std_normal()**2)
  if (a2 .gt. 0.d0) then
    write(*,*) i, a2, '# OK #'
  else
    write(*,*) i, a2, '# BUG #'
  end if 
end do
end program bug 

Output with 32 bit compiler:

 Vern : FTN95 v9.03.0
 Options : CHECKMATE;ERROR_NUMBERS;LINK;NO_BANNER;UNLIMITED_ERRORS;VS7;

 OK with intermediate value
            1         0.534344548168             0.285524096156
            2        -0.124214116238           1.542914667290E-02
            3         0.986227639562             0.972644957036
            4         -1.01574729998              1.03174257742
            5         0.876187316155             0.767704212991
            6         0.689662288046             0.475634071553
            7         0.236283819468           5.583004334239E-02
            8        -0.172203452071           2.965402890512E-02
            9          1.94609445045              3.78728361007
           10      -7.806535318409E-02         6.094199367757E-03

 Bug appears without the intermediate value
            1         0.942003631129    # OK #
            2        -0.647829889894    # BUG #
            3         0.247587151446    # OK #
            4         0.590547558638    # OK #
            5         0.337313028513    # OK #
            6         0.612344923211    # OK #
            7         0.233890118419    # OK #
            8        -0.526714233628    # BUG #
            9      -9.626477024171E-02  # BUG #
           10         0.296546433141    # OK #

Hi Ken,

look at this program !!

module random_numbers
implicit none
integer, parameter :: dp = kind(1.d0)
real(dp), parameter :: pi = 4.d0*atan(1.d0)
contains

real(dp) function random_std_uniform()
! Function to generate a random number from the standard uniform distribution.
 real(dp) :: x
 do
!z   call random_number(random_std_uniform)
   call random_number(x)
   if (x .gt. tiny(1.0)) exit       !  ?? tiny(1.0_dp)
 end do
 random_std_uniform = x
end function random_std_uniform

real(dp) function random_std_normal()
! Function to generate a random number from the standard normal distributon. 
 real(dp) :: x,y,r
 x = random_std_uniform()
 y = random_std_uniform()
 r = sqrt(-2.d0*log(x)) * cos(2.d0*pi*y)
 write (*,*) 'r =',r
 random_std_normal = r
end function random_std_normal

end module random_numbers

program bug
use random_numbers
use iso_fortran_env
implicit none
integer i
real(dp) a, a2
write(*,*) 'Vern : ',compiler_version ()
write(*,*) 'Options : ',compiler_options ()
write(*,*)
! This do loop runs as expected.  a can be positive or negative, a**2 is always positive
write(*,*) 'OK with intermediate value'
do i = 1, 10
  a = random_std_normal()
  a2 = (a**2)
  write(*,*) i, a, a2
end do
write(*,*)
! This do loop which does not have an intermediate set does not run as expected with the
! 32 bit compiler, it will on occasion return a negative number.
write(*,*) 'Bug appears without the intermediate value'
do i = 1, 10
  a2 = (random_std_normal()**2)
  if (a2 .gt. 0.d0) then
    write(*,*) i, a2, '# OK #'
  else
    write(*,*) i, a2, '# BUG #'
  end if
end do
end program bug

I ran this program with Win32 ChechMate It looks like FTN95 /32 replaces a2 = (random_std_normal()**2) by a2 = random_std_normal() * random_std_normal()

ie two calls to random_std_normal()

This must be an optimisation to ** but gives an invalid result ? not sure why ?

ftn95 /64 does not do this

Steping through sdbg appeared to show the problem.

To stop the double evaluation, I also tried : a2 = ( ( random_std_normal() ) **2 )

in the hope that it would treat '( fn )' as a computed value, but the compilation of the revised statement still called the fn twice.

Converting : a2 = ( ( random_std_normal() ) **2 )

appears to become : a2 = ( random_std_normal() * random_std_normal() )

which is wrong to me.

For ** 2, FTN95 should evaluate the function to a first register (reg_f) then copy the first register value to other register for multiply.

I don't see that double evaluation of function is a valid interpretation or an efficiency conversion from ** to reg_f * rer_f

FTN95 /64 only has 1 evaluation of the function, which appears to work for both loops.

John,

It took me a while to whittle the program that these module functions came from down to something that reproduced the error, and I never got as far as stepping through the reproducer one line at a time. Thanks for continuing the detective work!

The oddity was that chi-squared random variables could be negative! I only spotted this when comparing histograms of the variables produced by the 32 bit and 64 bit compilers. These are calculated with the following simple function:

real(dp) function chi_squared_rand(k)
! Function to generate a random number from the chi_squared distribution with K degrees of freedom
integer, intent(in) :: k
integer  :: j
  chi_squared_rand = 0.d0
  do j = 1, k, 1
    chi_squared_rand = chi_squared_rand + random_std_normal()**2
  end do
end function chi_squared_rand

As you can see it includes random_std_normal()**2

random_std_normal()**2 is definitely not the same as random_std_normal()*random_std_normal() and the 32 bit compiler should not call the function twice.

It does not occur with the 64 bit FTN95, or the other compilers I have access to.

I have being trying to make code as compact as possible, without unnecessary intermediate values. Got caught out of this occasion, and now I wondering where else I many have function()**2 in a 32 bit executable.

All part of the fun of coding I guess!

Ken

Here is a simplified reproducer.

program p
implicit none
real*8 :: a, a2
integer :: i

write(*,*) 'a =random_pmhalf(),  followed by a2 = a**2'

do i = 1, 10
  a = random_pmhalf()
  a2 = a**2
  if (a2 .ge. 0.d0) then 
    write(*,*) i, a2, 'PASS'
  else
    write(*,*) i, a2, 'FAIL'
  end if
end do   

write(*,*)
write(*,*) 'Failure with  random_pmhalf()**2'

do i = 1, 10
  a2 = random_pmhalf()**2
  if (a2 .ge. 0.d0) then
    write(*,*) i, a2, 'PASS'
  else
    write(*,*) i, a2, 'FAIL'
  end if
end do   

contains

real*8 function random_pmhalf()
! return a uniform random number in the range -0.5 to + 0.5
real*8 x
  call random_number(x)
  random_pmhalf = x - 0.5d0
end function  random_pmhalf

end program p

Which returns with the 32 bit compiler:

 a =random_pmhalf(),  followed by a2 = a**2
            1       2.995336023542E-02  PASS
            2         0.124266994006    PASS
            3       2.428571335641E-02  PASS
            4       5.590783884962E-02  PASS
            5       1.311532721418E-02  PASS
            6         0.244499685923    PASS
            7       1.775696581047E-03  PASS
            8       4.215554490142E-03  PASS
            9       5.499879256713E-02  PASS
           10         0.115877952355    PASS

 Failure with  random_pmhalf()**2
            1      -3.799733823335E-02  FAIL
            2         0.121547933359    PASS
            3       7.539772367957E-02  PASS
            4         0.170505953202    PASS
            5       8.266413839442E-02  PASS
            6       6.941677320227E-02  PASS
            7      -3.830583412939E-03  FAIL
            8      -6.198978447679E-04  FAIL
            9       4.229278825887E-03  PASS
           10      -1.743566815011E-02  FAIL

Many thanks (Ken) for the bug report and (John) for the investigation. This has now been fixed for the next release of FTN95.