 |
forums.silverfrost.com
Welcome to the Silverfrost forums
|
| View previous topic :: View next topic |
| Author |
Message |
Kenneth_Smith
Joined: 18 May 2012
Posts: 730
Location: Hamilton, Lanarkshire, Scotland.
|
 Posted: Sat Oct 12, 2024 5:39 pm Post subject: 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.
| Code: | 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:
| Code: | 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 # |
|
|
| Back to top |
|
 |
JohnCampbell
Joined: 16 Feb 2006
Posts: 2593
Location: Sydney
|
| Posted: Sun Oct 13, 2024 3:33 am Post subject: |
|
|
Hi Ken,
look at this program !!
| Code: | 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 |
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2593
Location: Sydney
|
| Posted: Sun Oct 13, 2024 9:05 am Post subject: |
|
|
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. |
|
| Back to top |
|
|
Kenneth_Smith
Joined: 18 May 2012
Posts: 730
Location: Hamilton, Lanarkshire, Scotland.
|
| Posted: Sun Oct 13, 2024 11:41 am Post subject: |
|
|
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:
| Code: | 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 |
|
| Back to top |
|
|
Kenneth_Smith
Joined: 18 May 2012
Posts: 730
Location: Hamilton, Lanarkshire, Scotland.
|
| Posted: Sun Oct 13, 2024 12:41 pm Post subject: |
|
|
Here is a simplified reproducer.
| Code: | 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:
| Code: | 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 |
|
|
| Back to top |
|
|
PaulLaidler
Site Admin
Joined: 21 Feb 2005
Posts: 8067
Location: Salford, UK
|
| Posted: Mon Oct 14, 2024 8:10 am Post subject: |
|
|
| Many thanks (Ken) for the bug report and (John) for the investigation. This has now been fixed for the next release of FTN95. |
|
| Back to top |
|
|
|
|
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|
FAQ
Search
Memberlist
Usergroups
Register
Profile
Log in to check your private messages
Log in
