Re: Scaling 64-bit random numbers



On 5/11/2011 9:46 PM, Gordon Sande wrote:
On 2011-05-11 21:49:37 -0300, Michael Trover said:

On 5/11/2011 7:13 PM, remek wrote:
Hello,

I have been looking for a good 64-bit pseudo-random number generator
for a long time and it looks like the best way to obtain something
satisfactory and reliable is to implement a generator oneself. For
this purpose, Marsaglia's KISS generator adapted for 64-bit numbers
seems to be an excellent choice:
http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/a85bf5f2a97f5a55


However,

this algorithm generates random 64-bit integers, whereas I
need double random numbers between 0 and 1. I am therefore wondering
how to scale these double integers appropriately.

A simple way would be to scale numerically by dividing the integer
generated by KISS() by 2*(2^63-1) and add 0.5 as follows:

real(8) function runif1()
implicit none
runif1 = 0.5_8 + KISS() * 0.542101086242d-19
end function runif1

Another approach would be to use the random bit sequence of the
integer generated by KISS() and to convert the 52 rightmost bits into
the mantissa of a double real number between 1 and 2, and finally
subtract 1 to obtain the correct range. This could be done by applying
logical masks as follows:

real(8) function runif2()
implicit none
integer(8) :: x
x = KISS()
x = iand(x, z'FFFFFFFFFFFFF')
x = ieor(x, z'3FF0000000000000')
runif2 = transfer(x, runif2) - 1.0_8
end function runif2

The first approach (runif1) troubles me a little because of numerical
accuracy. The second approach (runif2) seems to work fine on my
computer, but I am wondering if that makes sense and how portable it
is.

Any help on this would be greatly appreciated!

Thanks a lot in advance,
remek

The 'correct' way is as follows. Assuming your 64 bit generator
generates all whole numbers with 63 bits (the 64th bit is the sign
bit), and r is a double precision real number then
nummax = 2.0D0**63
r = KISS()/nummax.
Then r is between 0 and 1 and has a distribution defined by KISS().
This is the same as masking out the top 53 bits of the whole number
and applying it to r.

MikeT

Notice that the reals between 0.5 and 1.0 will be more widely spaced
than those
between 0.25 and 0.5 and so on. Once the integer has enough leading
zeroes so
that it will not loose low order bits this changing spacing will calm down.
Whether this is an important issue can only be answered by the original
poster.
One guesses the question arose because of a concern for rare events so the
question is not trivial.

This quantization issue is not the same as the question of whether the
random
number generator has an adequate period.


Hey, Gordon. I posted a subsequent reply before I got yours. You might want to look at that as my answer.

The quantization effect is a property of the finite representation of real numbers and is not unique or affected by the distribution of integers used to generate them. If the OP is concerned about this quantization effect then he is concerned about it everywhere real numbers are used in this real number representation.

MikeT
.



Relevant Pages

  • Scaling 64-bit random numbers
    ... I have been looking for a good 64-bit pseudo-random number generator ... Marsaglia's KISS generator adapted for 64-bit numbers ... realfunction runif1() ... realfunction runif2() ...
    (comp.lang.fortran)
  • Re: Scaling 64-bit random numbers
    ... On 5/11/2011 7:13 PM, remek wrote: ... I have been looking for a good 64-bit pseudo-random number generator ... Marsaglia's KISS generator adapted for 64-bit numbers ... realfunction runif2() ...
    (comp.lang.fortran)
  • Re: Scaling 64-bit random numbers
    ... I have been looking for a good 64-bit pseudo-random number generator ... Marsaglia's KISS generator adapted for 64-bit numbers ... realfunction runif1() ... realfunction runif2() ...
    (comp.lang.fortran)
  • Re: integer random number generator in given range
    ... Marsaglia's KISS RNGs typically return integers, ... you will get totally wrong results when the generator ... and double precision values. ...
    (comp.lang.fortran)