長周期乱数（L’Ecuyer）のPythonプログラム
合同式を用いる乱数発生法は、Pythonでは整数型に桁数の制限が文法上は無いので、実行速度を無視すれば、式のままプログラムを書くことができる。例えば、L’Ecuyer の方法（cf. Eubank & Kupresanin, 2012. Statistical computing in C++ and R.）は次の２つの式
に基づくものであるが、これはPythonでは次のように
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
簡単に書くことができる。具体的には、ファイルRNGen.pyにクラス型として宣言したが、上の説明と直接関係する部分は以下のようである。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ 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)
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
乱数生成のオブジェクトは、初期化時に６個の値を設定することができる。初期化時の
、
、
の値が全てが０ということではなく、かつ
より小さい、および初期化時の
、
、
の値が全てが０ということではなく、かつ
より小さいとき、乱数の周期は約
である。本ウェブサイトのプログラムファイルRNGen.pyなどは、圧縮ファイルrngprgs.zipにまとめたが、ダウンロード解凍して、利用者の責任において、自由に利用されたい。
乱数生成のクラスRNp2to191の使用例のプログラムを示す。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ 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()))
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のプログラムでは
from RNGen import RNp2to191 as rng
とクラスRNp2to191をimportしているので、
rn1 = rng()
rn2 = rng(10, 20, 30, 40, 50, 60)
と名前rngによってクラス型が示されている。オブジェクトrn1の初期値はデフォルト値であり、rn2は初期値が指定されている、このことは
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))
による出力によって確認することができる。
プログラム後半における関数jump(n)は、乱数生成において途中の
個の乱数を飛び越えるものである。これは、以下の関係を利用している。
例えば、式（１）の関係式は行列を使って以下のように表すことができる。
したがって、
とおけば、
と書けるので、乱数の生成は行列
を左から掛けることとして表すことができる。したがって、m回後の乱数は
を左から掛けることによって得られる。関数jump(n)は、
を求めて掛けている。これは、行列を２乗することをn回繰り返すことによって得られる。関数jumpのコードは以下のようである。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ 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
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のプログラムCheckRN.pyを実行すると図１に示すような結果になる。
図１
2つの乱数オブジェクトrn1とrn2の初期値が出力され、その後、関数jumpの効果を見る出力が表示されている。
関数jumpによって乱数列を分割して部分列に分け、それぞれを独立した乱数として別個のオブジェクトで表せば、独立した乱数列を複数個使うことができる。どれくらいの部分列に分ければよいかの目安を得るため、
個の乱数を生成する時間を測るプログラムCheckTime.pyを作成した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ 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))
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のプログラムの実行結果を図２に示す。
図２
乱数を個生成するのに掛かる日数が約
日であるとの結果である。現実問題として、個の乱数を使い切ることはない。関数jump(100)で部分列に区切って使えば、独立な乱数列が得られることになる。乱数全体の周期はであるので、個の乱数で区切った場合、
個の部分列が使用可能である。