A Random Number Generator Program Purely in Python
A random number generator with linear
congruential methods can be easily written purely in Python when the speed is
ignored, because, grammatically, the type integer in Python has no limits with
respect to its bit length. The generator by L’Eculyer
(cf. Eubank & Kupresanin, 2012, Statistical computing in C++ and R.) uses
the following congruential equations:
The rules of Equations (1) and (2) can be simply written in
Python as follows:
tx1
= (1403580 * self.x11 - 810728 * self.x12) % self.c209
self.x12 = self.x11
self.x11 = self.x10
self.x10 = tx1
tx2
= (527612 * self.x20 - 1370589 * self.x22) % self.c22853
self.x22 = self.x21
self.x21 = self.x20
self.x20 = tx2
In the program RNGen.py, the part, which directly correspond
to the above equations, is written as follows.
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ From the file
RNGen.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#
#
A progtram purely in Python
#
of
#
Pierre L'Ecuyer's algorithm of random number generation
#
#
Yasuharu Okamoto,2016.09
#
# cf.
R. L. Eubank and A. Kupresanin (2012). "Statistical computing in C++ and
R".
#
class RNp2to191:
c209 = 2**32 - 209
c22853
= 2**32 - 22853
def
__init__(self,*t):
if
(len(t) == 0):
self.x10 = 3
self.x11 = 2
self.x12 = 1
self.x20 = 6
self.x21 = 5
self.x22 = 4
else:
if (len(t) != 6):
print("len(t) = %d"%(len(t)))
raise "Number of parameters Error !"
else:
self.x10 = t[0]
self.x11 = t[1]
self.x12 = t[2]
self.x20 = t[3]
self.x21 = t[4]
self.x22 = t[5]
#
#
Uniform distribution, 0 < uni() < 1
#
def
uni(self):
tx1
= (1403580 * self.x11 - 810728 * self.x12) % RNp2to191.c209
self.x12 = self.x11
self.x11 = self.x10
self.x10 = tx1
tx2
= (527612 * self.x20 - 1370589 * self.x22) % RNp2to191.c22853
self.x22 = self.x21
self.x21 = self.x20
self.x20 = tx2
z =
(self.x10 - self.x20) % RNp2to191.c209
if (z > 0):
return z / (RNp2to191.c209 + 1.0)
else:
return RNp2t0191.c209 / (RNp2to191.c209 + 1.0)
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
When constructing an object of the above class RNp2to191,
initial values of the seeds can be set. When the initial values satisfy the
following conditions
, 
and 
are not all zero, and they are less than 
and
, 
and 
are not all zero, and they are less than 
,
The period of the generator is about 
. This period
is long enough so that independent sub-chains can be easily obtained, which is
explained in the latter half of the website.
Programs, which are explained in this website, are archived
into the archived file rngprgs.zip, and the programs
in the archived file rngprgs.zip can be used freely by
the user on his/her responsibility, although all rights are reserved by the
author.
An example of using the random number generator class
RNp2to191 is shown in the following file CheckRN.py.
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ CheckRN.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
from RNGen import RNp2to191 as rng
rn1 = rng()
rn2 = rng(10, 20, 30, 40, 50, 60) # The seeds
are set
print("rn1 = %d, %d, %d, %d, %d, %d"%(rn1.x10,
rn1.x11, rn1.x12,rn1.x20,rn1.x21,rn1.x22))
print("rn2 = %d, %d, %d, %d, %d, %d"%(rn2.x10,
rn2.x11, rn2.x12,rn2.x20,rn2.x21,rn2.x22))
rn_0 = rng()
# No jumps
rn_2 = rng()
rn_2.jump(2)
# 2**2 jumps
rn_22 = rng()
rn_22.jump(2) # 2**2 jumps
rn_22.jump(2) # 2**2 jumps, again
rn_4 = rng()
rn_4.jump(4)
# 2**4 jumps
for i in range(0, 20):
print("%d %20.15f %20.15f %20.15f %20.15f"%(i, rn_0.uni(),
rn_2.uni(), rn_22.uni(), rn_4.uni()))
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
In the above program, the class RNp2to191 is imported as
follows:
from RNGen import RNp2to191 as rng
Since the name of the class RNp2to191 is denoted as rng, we
can create objects of the class RNp2to191 by the following code:
rn1 = rng()
rn2 = rng(10, 20, 30, 40, 50, 60)
The initial values of the object rn1 are default values, and
those of the object rn2 are set by the parameter values. Those values set as
initial ones can be checked by the following codes.
print("rn1 = %d, %d, %d, %d,
%d, %d"%(rn1.x10, rn1.x11, rn1.x12,rn1.x20,rn1.x21,rn1.x22))
print("rn2 = %d, %d, %d, %d,
%d, %d"%(rn2.x10, rn2.x11, rn2.x12,rn2.x20,rn2.x21,rn2.x22))
Calling the function jump(n) causes a jump over 
random numbers in the series to be
generated. This uses the following mathematical relation.
Equation (1) can be written by matrix notation as follows
Set
then we have
Hence, the generation of the random number can be
represented by multiplication of the matrix 
.
Multiplication of 
generates the random number after m
steps. The function jump(n) calculates 
, which is
obtained by squaring of matrices n times, then jump over random numbers to be generated. The
source code of the function jump is as follows.
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ From the file
RNGen.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#
#
Jump over 2**n numbers to be generated
#
def
jump(self, n):
if
(n > 0):
A = [[0, 1403580, -810728],
[1, 0, 0],
[0, 1, 0]]
B = [[527612, 0, -1370589],
[1, 0, 0],
[0, 1, 0]]
C = [[0, 0, 0],
[0, 0, 0],
[0, 0, 0]]
for i in range(0, n):
for
j in range(0, 3):
for k in range(0, 3):
tv = 0
for h in range(0, 3):
tv = tv + A[j][h] * A[h][k]
C[j][k] = tv
for j in
range(0, 3):
for k in range(0, 3):
A[j][k] = C[j][k] % RNp2to191.c209
for j in range(0, 3):
for k in range(0, 3):
tv = 0
for h in range(0, 3):
tv = tv + B[j][h] * B[h][k]
C[j][k] = tv
for j in range(0, 3):
for k in range(0, 3):
B[j][k] = C[j][k] % RNp2to191.c22853
t0 = A[0][0] * self.x10 + A[0][1] * self.x11 + A[0][2] * self.x12
t1 = A[1][0] * self.x10 + A[1][1] * self.x11 + A[1][2] * self.x12
t2 = A[2][0] * self.x10 + A[2][1] * self.x11 + A[2][2] * self.x12
self.x10 = t0 % RNp2to191.c209
self.x11 = t1 % RNp2to191.c209
self.x12 = t2 % RNp2to191.c209
t0 = B[0][0] * self.x20 + B[0][1] * self.x21 + B[0][2] * self.x22
t1 = B[1][0] * self.x20 + B[1][1] * self.x21 + B[1][2] * self.x22
t2 = B[2][0] * self.x20 + B[2][1] * self.x21 + B[2][2] * self.x22
self.x20 = t0 % RNp2to191.c22853
self.x21 = t1 % RNp2to191.c22853
self.x22 = t2 % RNp2to191.c22853
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
The results of execution of the program CHeckRN.py is shown
in Figure 1.
Figure 1
First, the initial values of the objects rn1 and rn2 are
shown, then, effects of calling the function jump are displayed. Sequences of
random numbers are adequately shifted.
Construction
of Independent Chains
We can use multiple independent sequences of random numbers,
when the original sequence is divided into subsequences, which can be obtained
by the function jump. What length of subsequences is required. I uses 
.
The following program CheckTime.py checks the time, which
will be took by generation of random numbers.
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ CheckTime.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import time
from RNGen import RNp2to191 as rng
print("I am working...\n")
rn = rng()
t0 = time.time()
n = 2**23
for i in range(0,n):
rn.uni()
t1 = time.time()
print("Generation of 2**23 random numbers took about %dsec."%(t1
- t0))
t = (t1 - t0) * ((2.0**100) / (2.0**23))
days = t / (24.0 * 60.0 * 60.0)
print("Generation of 2**100 random numbers will take
about %dsec.,"%t +
"i.e.,\n
about%8.1edays."%days)
print("2**23 = %d"%(2**23))
print("2**100 = %d"%(2**100))
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
The results are shown in Figure 2.
Figure 2
It took about 
days to generate random numbers. In reality, random numbers cannot be exhausted.
Partitioning of the sequence by calling jump(100), if necessary several times
according to the required number of subsequences, provides independent
subsequences. Since 
, we can use
at most 
independent subsequences.
Yasuharu Okamoto