jeudi 16 juin 2016

Upper bound of random number generator

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