前回のコードは、元にしたプログラムがゲーム用に最適化されちゃっていて、もとの論文と比較しにくかったので、乱数生成部分を書き直しました。
WELL512a と WELL1024a, はC版と同じ結果が返ることを確認しています。
それ以外は、たぶん動くけど未検証なので、コメントアウトしてあります。
（誰か検証してください！）
あと、jumpahead()は、未実装にしました。seed()使って下さい。
本来の意味だと、数学的に離れた位置にジャンプするべきなのですが、
やり方がわからないんで、前はseed()を使ってお茶を濁してました。
これは、かえって混乱のもとなのでやめました。
# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# WELL(512a, 1024a)のPython実装
#
# 元の論文との対応付けのし易さを最重視して実装しています。
# 状態の初期化にSHA512を使うのがオリジナルのアイディアです。
# 512より多い状態を初期化する際を考慮して、
# SHA512の繰り返しとXorShiftのXORで状態を初期化しています。
#
# パフォーマンスを求める方は、CやC++版を使って下さい。
#
# このコードを作成するにあたり、参考にしたのは、以下です
#
# - http://d.hatena.ne.jp/fgshun/20090724/1248446400
# - http://www.iro.umontreal.ca/~panneton/WELLRNG.html
# - http://www3.ocn.ne.jp/~harase/megenerators.html
# - http://bitbucket.org/sergiu/random/src/tip/well.hpp
# - http://en.wikipedia.org/wiki/Well_Equidistributed_Long-period_Linear
#
# このコードのライセンスは、MITでお願いします。
#
# ---------------------------------------------------------------------------
from random import Random as _Random
from os import urandom as _urandom
from binascii import hexlify as _hexlify
from hashlib import sha512 as _sha512
from itertools import cycle as _cycle, izip as _izip
# ---------------------------------------------------------------------------
# 各WELLのベースクラス定義
#
# random_int32()が未実装なので、単独では使用できません。
# init_subclassと_M1～_M6を使用して、random_int32()を初期化して使って下さい。
# ---------------------------------------------------------------------------
def _xor128(x=123456789, y=362436069, z=521288629, w=88675123):
while True:
t = x^((x << 11) & 0xffffffff)
x = y
y = z
z = w
w = (w ^ (w >> 19)) ^ (t ^ (t >> 8))
yield w
class _WELLBase(_Random):
_W = 32
_MASK = 2 ** _W - 1
def __init__(self, seed=None):
self.state = self._gen_seed(seed)
self._check_zerostate()
self.gauss_next = None
self.index = 0
def _check_zerostate(self):
if self.state == [0] * self._R:
self.state = self._gen_seed(self.VERSION)
@classmethod
def _gen_seed(cls, a=None):
if a is None:
try:
a = int(_hexlify(_urandom(16)), 16)
except NotImplementedError:
import time
a = int(time.time() * 256)
if not isinstance(a, basestring):
a = str(a)
hashed_states = _sha512(a).hexdigest()
states = [int(hashed_states[i: i+8], 16) for i in range(0, 128, 8)]
return [xor128 ^ sha512 for (i, xor128, sha512) in _izip(range(cls._R), _xor128(), _cycle(states))]
def seed(self, a=None):
self.__init__(a)
def random_int32(self):
raise NotImplementedError
def random(self):
return self.random_int32() * 2.3283064365386963e-10
def random53(self):
return self.getrandbits(53) * (1.0 / 9007199254740992.0)
def getstate(self):
return (self.VERSION, tuple(self.state),
self.gauss_next, self.index)
def setstate(self, state):
version = state[0]
if version == self.VERSION:
self.state = list(state[1])
self.gauss_next = state[2]
self.index = state[3]
else:
raise ValueError("state with version %s passed to "
"setstate() of version %s" %
(version, self.VERSION))
def jumpahead(self, n):
raise NotImplementedError
def getrandbits(self, k=32):
ret = self.random_int32()
for _ in range((k - 1) // 32):
ret = (ret << 32) ^ self.random_int32()
return ret >> ((- k) % 32)
@classmethod
def init_random_int32(cls, T0, T1, T2, T3, T4, T5, T6, T7, m1, m2, m3, me_mask=0):
MASKU = cls._MASK >> (cls._W - cls._P)
MASKL = cls._MASK << cls._P
r = cls._R
def random_int32(self):
s = self.state
newV1 = i = self.index
newV0 = i_VRm1 = (i - 1) % r
i_VRm2 = (i - 2) % r
vm2 = s[(i + m2) % r]
z0 = (s[i_VRm1] & MASKL) | (s[i_VRm2] & MASKU);
z1 = T0(s[i]) ^ T1(s[(i + m1) % r])
z2 = T2(vm2 ) ^ T3(s[(i + m3) % r])
s[newV1] = z3 = z1 ^ z2
s[newV0] = z4 = T4(z0) ^ T5(z1) ^ T6(z2) ^ T7(z3)
self.index = i_VRm1
return z4 ^ (vm2 & me_mask)
cls.random_int32 = random_int32
@classmethod
def _output_after_tempering(cls, b, c):
random_int32 = cls.random_int32
def random_int32_after_tempering(self):
y = random_int32(self)
y = y ^ ((y << 7) & b)
y = y ^ ((y << 15) & c)
return y
cls.random_int32 = random_int32_after_tempering
@classmethod
def _shift(cls, k):
MASK = cls._MASK
if k >= 0:
return lambda v: (v >> k)
else:
k_abs = abs(k)
return lambda v: (v << k_abs) & MASK
@classmethod
def _M0(cls):
return lambda v: 0
@classmethod
def _M1(cls):
return lambda v: v
@classmethod
def _M2(cls, t):
return cls._shift(t)
@classmethod
def _M3(cls, t):
t_shift = cls._shift(t)
return lambda v: v ^ t_shift(v)
@classmethod
def _M4(cls, a):
w1 = cls._W - 1
return lambda x: (x >> 1) ^ ((x & 1) * a)
@classmethod
def _M5(cls, t, b):
t_shift = cls._shift(t)
return lambda v: v ^ (t_shift(v) & b)
@classmethod
def _M6(cls, r, a, ds, dt):
q_shift = cls._shift( - r)
wq_shift = cls._shift(cls._W - r)
return lambda x: ((q_shift(x) ^ wq_shift(x)) & ds) ^ (a if (x & dt) != 0 else 0)
# ------------------------------------------------------
# 個別のWELLの定義
# ------------------------------------------------------
class WELL512a(_WELLBase):
VERSION = u"WELL512a ver0"
_K = 512 # not use
_R = 16
_P = 0
WELL512a.init_random_int32(
m1=13, m2=9, m3=5,
T0=WELL512a._M3(-16), T4=WELL512a._M3(-2),
T1=WELL512a._M3(-15), T5=WELL512a._M3(-18),
T2=WELL512a._M3(11), T6=WELL512a._M2(-28),
T3=WELL512a._M0(), T7=WELL512a._M5(-5, 0xda442d24))
class WELL1024a(_WELLBase):
VERSION = u"WELL1024a ver0"
_K = 1024 # not use
_R = 32
_P = 0
WELL1024a.init_random_int32(
m1=3, m2=24, m3=10,
T0=WELL1024a._M1(), T4=WELL1024a._M3(-11),
T1=WELL1024a._M3(8), T5=WELL1024a._M3(-7),
T2=WELL1024a._M3(-19), T6=WELL1024a._M3(-13),
T3=WELL1024a._M3(-14), T7=WELL1024a._M0())
# 以下、未検証のためコメントアウト
#
#class WELL1024b(_WELLBase):
# VERSION = u"WELL1024b ver0"
# _K = 1024 # not use
# _R = 32
# _P = 0
#
#WELL1024b.init_random_int32(
# m1=22, m2=25, m3=26,
# T0=WELL1024b._M3(-21), T4=WELL1024b._M3(-14),
# T1=WELL1024b._M3(17), T5=WELL1024b._M3(-21),
# T2=WELL1024b._M4(0x8bdcb91e), T6=WELL1024b._M1(),
# T3=WELL1024b._M3(15), T7=WELL1024b._M0())
#
#class WELL19937a(_WELLBase):
# VERSION = u"WELL19937a ver0"
# _K = 19937 # not use
# _R = 624
# _P = 31
#
#WELL19937a.init_random_int32(
# m1=70, m2=170, m3=449,
# T0=WELL19937a._M3(-25), T4=WELL19937a._M1(),
# T1=WELL19937a._M3(27), T5=WELL19937a._M3(-9),
# T2=WELL19937a._M2(9), T6=WELL19937a._M3(-21),
# T3=WELL19937a._M3(1), T7=WELL19937a._M3(21))
#
#class WELL19937b(_WELLBase):
# VERSION = u"WELL19937b ver0"
# _K = 19937 # not use
# _R = 624
# _P = 31
#
#WELL19937b.init_random_int32(
# m1=203, m2=613, m3=123,
# T0=WELL19937b._M3(7), T4=WELL19937b._M3(-19),
# T1=WELL19937b._M1(), T5=WELL19937b._M2(-11),
# T2=WELL19937b._M3(12), T6=WELL19937b._M3(4),
# T3=WELL19937b._M3(-10), T7=WELL19937b._M3(-10))
#
#class WELL19937c(WELL19937a):
# VERSION = u"WELL19937c ver0"
#
#WELL19937c._output_after_tempering(0xe46e1700, 0x9b868000)
#
#class WELL44497a(_WELLBase):
# VERSION = u"WELL44497a ver0"
# _K = 44497 # not use
# _R = 1391
# _P = 15
#
#WELL44497a.init_random_int32(
# m1=23, m2=481, m3=229,
# T0=WELL44497a._M3(-24), T4=WELL44497a._M1(),
# T1=WELL44497a._M3(30), T5=WELL44497a._M3(20),
# T2=WELL44497a._M3(-10), T6=WELL44497a._M6(9, 0xb729fcec,
# 0xfbffffff, 0x00020000),
# T3=WELL44497a._M2(-26), T7=WELL44497a._M1())
#
#class WELL44497b(WELL44497a):
# VERSION = u"WELL44497b ver0"
#
#WELL44497b._output_after_tempering(0x93dd1400, 0xfa118000)
#
# EOF
イッシーの利用目的だと、WELL512aでおそらく十分ですので、
要望や助力がない限りは、これ以上の実装例は増やさない予定です。
英語読めなくても、論文とソースを見比べていけば、なんとななるはずです。
（イッシーがそうです）
あと、init_random_int32()の、me_maskは、高次最適化の為のオプションです。
これはそのうち対応したいなぁ。
0 件のコメント:
コメントを投稿