Re: double precision RNGs
George Marsaglia
To produce, in the unit interval 0<=x<1, a double
precision random variable x under the IEEE 754
standard, you need to provide, in effect, the
floating point version of a ratio k/2^53
for some random k of 53 or fewer bits.
It is possible, as in the version I reported here, (Nov 18)
based on ``The 64-bit universal RNG", (2004)
Statistics and Probability Letters v.66 No. 2, 183-187,
to do this using only floating point addition
and subtraction so that exactly the same results
are produced for a variety of platforms and
programming languages.
But if you have a good routine for producing
random integers, it is relatively easy to
use them to produce IEEE double precision x's:
For example, if rand() produces unsigned 32-bit
integers, then
x=(2097152.*rand( )+(rand( )>>11))/9007199254740992.;
will work in C. In Fortran, which lacks the
ability to handle unsigned 32-bit integers, any
31-bit RNG, say h(), may be used to provide
x=(4194304.*h( )+ishift(h( ),9))/9007199254740992.
(The floating point constants above are 2.^21, 2.^53
and 2.^22. Be sure to include their decimal points.)
Of course, if you want a good x from this method,
you need a good source of 31- or 32-bit random
integers. The rand( ) functions available with
most C compilers are not very good.
My favorite, unless I want to use one of the MWC
or CMWC RNGs with massive periods, but requiring
hundreds to thousands of seeds, is KISS( ), the
Keep It Simple Stupid RNG that combines congruential,
xorshift and lag-1 multiply-with-carry sequences:
unsigned long KISS( )
{static unsigned long
x=123456789,y=362436000,z=521288629,c=7654321;
unsigned long long t, a=698769069LL;
x=69069*x+12345;
y^=(y<<13); y^=(y>>17); y^=(y<<5);
t=a*z+c; c=(t>>32);
return x+y+(z=t);
}
Tke KISS RNG requires four seeds, the initial x,y,z and c.
It has period > 2^124 and passes all the tests I have
put to it, particularly the newer version of the
Diehard battery with some tough new tests,
http://www.csis.hku.hk/~diehard/
Since Fortran does not provide easy means to access
the top and bottom 32-bits of a 64-bit product,
a Fortran version uses two 16-bit MWC components:
function kiss( )
integer x,y,z,w
data x,y,z,w/123456789,362436069,521288629,916191069/
m(k,n)=ieor(k,ishft(k,n))
x=69069*x+1327217885
y=m(m(m(y,13),-17),5)
z=18000*iand(z,65535)+ishft(z,-16)
w=30903*iand(w,65535)+ishft(w,-16)
kiss=x+y+ishft(z,16)+w
return
end
To use the Fortran kiss( ) to produce the
IEEE double precision x, you could use
31 bits from kiss( ) to avoid getting negative
values in forming the 53 bits of the numerator
of x.
George Marsaglia
jacob navia
Here is KISS and frand for x86 assembly :-)
It makes more than 200 million frands per second
.data
_x:
.long 0x75bcd15
_y:
.long 0x159a55a0
_z:
.long 0x1f123bb5
_c:
.long 0x74cbb1
.text
_KISS:
movl $0x10dcd,%eax
mull _x
addl $12345,%eax
movl %eax,_x
movl _y,%edx
movl %edx,%ecx
shll $13,%ecx
xorl %ecx,%edx
movl %edx,%ecx
shrl $17,%ecx
xorl %ecx,%edx
movl %edx,%ecx
shll $5,%ecx
xorl %ecx,%edx
movl %edx,_y
movl $0x29a65ead,%eax
mull _z
addl _c,%eax
adc $0,%edx
movl %edx,_c
movl %eax,_z
addl _x,%eax
addl _y,%eax
ret
.align 2
_frand:
call _KISS
pushl $0
pushl %eax
call _KISS
fldl _$12
fildq (%esp)
shrl $11,%eax
fmulp %st,%st(1)
movl %eax,(%esp)
fildq (%esp)
faddp %st,%st(1)
addl $8,%esp
fdivl _$13
ret
.globl _frand
.data
.align 4
_$13:
.long 0x0,0x43400000
_$12:
.long 0x0,0x41400000