Prev: PR1ME OPERATING SYSTEM
Next: Pointer or (allocatable) array list ?
From: geo on 1 Jul 2010 05:27
I have found an interesting new
Random Number Generator (RNG) based on the prime

p = 2^137210-2^133114+1.

One of its most interesting features is that establishing
the period requires finding all of the prime factors of p-1,
and each of these famous giants of prime number theory:
Euler, Cunningham, Morrison, Brillhart,
Brent, Pollard, Lenstra, Selfridge
has discovered one or more of the prime factors of p-1:
p-1 = 2^133114*F_2*F_3*F_4*F_5*F_6*F_7*F_8*F_9*F_10*F_11,
with F_n=2^(2^n)+1 the nth Fermat number. See, for example,
http://www.prothsearch.net/fermat.html

If j is a random integer in 0<j<p then the binary expansion
of j/p will have period (p-1)/2 and its bits can well serve
as a source of random bits.

The bits in such an expansion may be generated---in reverse
order, either 32-bits at a time, or 64-bits at a time---by a
CSWB (Complimentary-Subtract-With-Borrow) procedure:

Given x_0,x_1,...,x_{4287},(32-bit seeds) and previously generated
x_{4288},x_{4289},...,x_{n-1} for n>4287, and a current `boro` of
0 or 1, these simple operations may be used to get
the next boro and the next random number x_n:
t=x_{n-4288}; h=x_{n-4160}+boro;
if t<h then boro=1 else boro=0;
return x_n=(h-t-1 mod b=2^32);
and keep boro for the next x_n.

The 32-bit version is based on the representation
p=b^4288-b^4160+1 with b=2^32, but setting B=2^64 and
p=B^2144-B^2080+1 leads to 64-bit integers by means of
the same instructions, except that, with 2144 64-bit seeds,
t=x_{n-2144}; h=x_{n-2080}+boro;

Note that x_n=h-t-1 will be automatically mod b=2^32
for most 32-bit CPUs, or mod B=2^64 for 64-bit CPUs,
whether for signed or unsigned integers.

Determining the new boro is particularly simple
for signed integers. In C: boro=(t<h);
For languages that permit only signed integers, it
is not as simple. Let s(x) be the sign bit of x,
easily obtained, for example, by a right shift:
if s(t)=s(h) then boro=s(t-h); else boro=s(h);

Interested readers are invited to provide C,C++,Fortran or
other implementations. All you need is an array, say X[4288]
for 32-bit or X[2144] for 64-bit CPUs, and return the next
unused element of the array, except that if none are left,
refill the array with the rules above, reset the array index
and return the first element. Suitable for signed or
unsigned, 32- or 64-bit integers, but signed may require the
more complicated method to determine the 'boro'.

The arrays X[4288] or X[2144] must be initially seeded with
32-bit or 64-bit integers and an initial boro seed of 0 or 1,
Having so many seed bits is desirable in applications such
as in Law or Gaming, where a minimal, often extremely large,
number of possible results are required, but a nuisance in most
applications where only a few dozen seed bits may be enough.
Most implementations would ask for a few seeds from which the
X[4288] or X[2144] arrays are filled by combining simple RNGs.

The CSWB sequence x_n=x_{n-4288}-x_{n-4160+boro mod 2^32
might do well on the Diehard tests, but I tend to favor
a KISS (Keep It Simple, Stupid) approach, combining it, via
addition mod b=2^32 or B=2^64, Congruential and Xorshift
sequences. Unlike most combinations, such a one would not
increase the period, as 2^32, 2^32-1, 2^64 and 2^64-1 all
divide p-1. But a period exceeding 10^3104 should serve.

George Marsaglia


 | 
Pages: 1
Prev: PR1ME OPERATING SYSTEM
Next: Pointer or (allocatable) array list ?