Hi there,
The intrinsic RANDOM_NUMBER in gfortran is said to be an
implementation of Marsaglia's KISS PRNG.
e.g. this post: http://gcc.gnu.org/ml/fortran/2007-07/msg00454.html
Looking at the implementation of random.c in libgfortran, the relevant
code is:
/* kiss_random_kernel() returns an integer value in the range of
(0, GFC_UINTEGER_4_HUGE]. The distribution of pseudorandom numbers
should be uniform. */
static GFC_UINTEGER_4
kiss_random_kernel(GFC_UINTEGER_4 * seed)
{
GFC_UINTEGER_4 kiss;
seed[0] = 69069 * seed[0] + 1327217885;
seed[1] = GFC_SL(GFC_SR(GFC_SL(seed[1],13),17),5);
seed[2] = 18000 * (seed[2] & 65535) + (seed[2] >> 16);
seed[3] = 30903 * (seed[3] & 65535) + (seed[3] >> 16);
kiss = seed[0] + seed[1] + (seed[2] << 16) + seed[3];
return kiss;
}
This coincides with the code found in http://www.fortran.com/kiss.f90
(of source unknown):
FUNCTION kiss ()
integer :: kiss
! The KISS (Keep It Simple Stupid) random number generator. Combines:
! (1) The congruential generator x(n)=69069*x(n-1)+1327217885, period
2^32.
! (2) A 3-****ft ****ft-register generator, period 2^32-1,
! (3) Two 16-bit multiply-with-carry generators, period
597273182964842497>2^59
! Overall period>2^123; Default seeds x,y,z,w.
! Set your own seeds with statement i=kisset(ix,iy,iz,iw).
!
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
contains
function m(k, n)
integer :: m, k, n
m = ieor (k, ishft (k, n) )
end function m
END FUNCTION kiss
However, this version does not coincide with Marsaglia's original
implementation given in the post
http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@[EMAIL
PROTECTED]
is
#include <stdio.h>
#define znew (z=36969*(z&65535)+(z>>16))
#define wnew (w=18000*(w&65535)+(w>>16))
#define MWC ((znew<<16)+wnew )
#define SHR3 (jsr^=(jsr<<17), jsr^=(jsr>>13), jsr^=(jsr<<5))
#define CONG (jcong=69069*jcong+1234567)
#define KISS ((MWC^CONG)+SHR3)
Notice the difference in the parameters for the multiple-with-carry
RNG MWC (36969 and 18000 in Marsaglia versus 18000 and 30903 in
gfortran, also in the wrong order), and in the CONG (69069 and 1234567
in Marsaglia versus 69069 and 1327217885 in gfortran). Anybody know
the story behind these changes? Anything published? I did an online
search but didn't find any reason. L'Ecuyer has studied Marsaglia's
original implementation and found it good, but I haven't seen any test
results on the gfortran version.
Final question, how do I recover the original integer sequence from
KISS, given the output from RANDOM_NUMBER. I don't think you just
divide the output by 2**32 do you?
Thanks for your help
Manuel
|
