I have recently posed questions to this group:
Form t = a*x+c in 64 bits, get top and bottom 32 bits of t
Form k = i + j and test for overflow
Such problems arise in certain kinds of random number
generation---multiply-with-carry and add-with-carry.
http://tbf.coe.wayne.edu:16080/jmasm/vol2_no1.pdf
I was hoping to develop RNGs that would produce
exactly the same sequence of 32-bit random integers
in Fortran as in C, but using only simple code.
It seems that it is not yet feasible to carry out,
in most Fortrans, using 32-bit integers,
the t=a*x+c operation for multiply-with-carry,
or the k=i+j overflow for add-with-carry.
However, add-with-carry is easy for 31-bit
integers.
I have therefore developed, and post here,
a version of my
KISS (Keep It Simple, Stupid) RNG,
which uses only
add, shift, exclusive-or and 'and' operations
to produces exactly the same 32-bit integer
output, which C views as unsigned and
Fortran views as signed integers,
each made up of exactly the same pattern
of 32 0's and 1's.
Furthermore, this common output sequence
has period greater that 10**36 and seems
to pass all the Diehard battery of tests
http://www.cs.hku.hk/~diehard/
and the tuftests in
http://www.jstatsoft.org/v07/i03/tuftests.pdf
/* C version */
static unsigned long t,x=123456789,y=362436069,z=21288629,w=14921776,c=0;
unsigned long KISS(){
x+=545925293;
y^=(y<<13); y^=(y>>17); y^=(y<<5);
t=z+w+c; z=w; c=(t>>31); w=t&2147483647;
return(x+y+w); }
! Fortran version, (using several statements/line via semicolons.)
function KISS()
integer x,y,z,w,c,t
data x,y,z,w,c/123456789,362436069,21288629,14921776,0/
m(y,k)=ieor(y,ishft(y,k))
x=x+545925293
y=m(m(m(y,13),-17),5)
t=z+w+c; z=w; c=ishft(t,-31); w=iand(t,2147483647)
KISS=x+y+w
return
entry kisset(ix,iy,iz,iw,ic)
x=ix; y=iy; z=iz; w=iw; c=ic
kisset=1
return
end
The output of both KISS() functions will be
a sequence of 32-bit integers with period greater
than 2**121 or 10**36:
period=576384491062058838*4294967296*4294967295
The value returned from KISS() can be any choice
of x op y op w, with op either '+' or 'exclusive-or'.
Each of the four possible choices produced good
results from the Diehard and tuftests, for a
variety of seed sets.
When b=2^31, b^2+b-1 is not prime, but factors into
610092078393289*7559
so the full period from KISS() requires that
for the seeds x,y,z,w,c:
x can be any 32-bit integer,
y can be any 32-bit integer not 0,
z and w any 31-bit integers not multiples of 7559
c can be 0 or 1.
I was surprised that three such simple
components, individually likely to fail
various tests of randomness, could
collectively, via x op y op w, produce
sequences that do as well in stringent tests
of randomness as do any I have previously
encountered.
That, with their long periods, simple
operations and resulting speed,
(84 million per second on my PC),
makes them worth considering.
Furthermore, by applying that last component,
add-with-carry, to a suitable pair of elements
from an array of 31-bit seeds z[1],z[2],...,z[n],
rather than just to the current z,w pair,
you can get periods on the order of 2**(31*n+64),
roughly the number of possible ways to initialize
x,y,c and the n 31-bit z's,
and have the same Fortran-C compatibility,
simplicity and virtually the same speed.
You are welcome to use one or both of the above,
and perhaps try assessing the randomness yourself.
(If you generate 10000 integers with
either the C or Fortran version listed
above, the last four values should be
Fortran C
199275006 199275006
86473693 86473693
-2085369775 2209597521
1298124039 1298124039
George Marsaglia