This is actually a follow up question of a previous one: Rounding of double precision to single precision: Forcing an upper bound
After what I thought was the solution of my problems with the answer of previous question, I tried running my program again and found that I had the same problem.
The Mersenne Twister implementation I'm using generates a signed 32 bits random integer. The guy who implemented the RNG made this function to generate a random double precision float in the range [0,1):
function genrand_real2()
double precision genrand_real2,r
integer genrand_int32
r=dble(genrand_int32())
if(r.lt.0.d0)r=r+2.d0**32
genrand_real2=r/4294967296.d0
return
end
And it works flawlessly, so following the suggestion in the previous question I used the following function to generate a random single precision float, in the range I thought would be [0,1):
function genrand_real()
real genrand_real, r
integer genrand_int32
r = real(genrand_int32())
if (r .lt. 0.0) r = r + 2.0**32
genrand_real = r / 4294967296.0
return
end
However I got the same error I got before, caused by a 1.0 number. So I wrote a small program to show that my genrand_real actually generates a 1.0, and found that I was right, and the 1.0 is generated. This causes the way I use to generate an integer in the range [1,MAX] (in this example [1,5]) to fail generating a value MAX+1.
i = 0
do while (.true.)
r = genrand_real()
if (r .gt. 0.99999) then
i = i + 1
print *, 'number is:', r
print *, 'conversion is: ', int(5*r)+1
endif
if (i .gt. tot_large) exit
enddo
My question is, why does it work for the double precision but not for the single precision float? I don't see a reason for it to fail since 2**32 fits in a single precision float. Also, what should I do to fix it? I thought about dividing the number by 2.0**32+1 instead of 2.0**32, but I'm not sure it's theoretically correct and that the numbers would be uniform.
Aucun commentaire:
Enregistrer un commentaire