WELL512aをPythonで実装してみる

擬似乱数生成法「WELL」の512ビット版をPythonで実装してみました！

このアルゴリズムは、まだマイナーのようですが、
有名なメルセンヌ・ツイスタ（Pythonの標準の乱数生成アルゴリズムでもあります）より新しく、より乱数の性質が良いとされているアルゴリズムです。

WELLについては、以下のサイトのメモがまとまっています。
http://lt140.blogspot.jp/2010/05/well.html
http://en.wikipedia.org/wiki/Well_equidistributed_long-period_linear

何気にまだちゃんと動作検証できていないのですが、勢い（だけ）でソースをのっけます。
jumpheadの動きが適当すぎるんで、もうちょっと変えます。

# -*- coding: utf-8  -*-
# WELL512aのPython実装
#
# 初期化にSHA512を使うのが特徴
#
# 参考にしたのは、以下です
#   - http://www.lomont.org/Math/Papers/2008/Lomont_PRNG_2008.pdf
#   - http://d.hatena.ne.jp/fgshun/20090724/1248446400
#
# このコードのライセンスは、MITでお願いします。
#

import random
from os import urandom as _urandom
from binascii import hexlify as _hexlify
from hashlib import sha512

_ALL_ZERO_LIST = [0] * 16


class WELL512(random.Random):

    VERSION = u"WELL512 ver0"

    def __init__(self, seed=None):
        _state = self._gen_seed(seed)

        if _state == _ALL_ZERO_LIST:
            self._state = self._gen_seed(self.VERSION)
        else:
            self._state = _state

        self.gauss_next = None
        self._index = 0

    @classmethod
    def _gen_seed(cls, a=None):
        if a is None:
            try:
                a = long(_hexlify(_urandom(16)), 16)
            except NotImplementedError:
                import time
                a = long(time.time() * 256) # use fractional seconds

        if not isinstance(a, basestring):
            a = str(a)

        seedhash = sha512(a)
        hashed_states = seedhash.hexdigest()
        return [long(hashed_states[i: i+8], 16) for i in range(0, 128, 8)]

    def seed(self, a=None):
        self.__init__(a)

    def _genrand_int32(self):
        index = self._index
        state = self._state

        a = state[index]
        c = state[(index + 13) & 15]
        b = (a ^ c ^ (a << 16) ^ (c << 15)) & 0xffffffff
        c = state[(index + 9) & 15]
        c ^= c >> 11
        a = state[index] = b ^ c
        d = a ^ ((a << 5) & 3661901092L)
        index = (index + 15) & 15
        a = state[index]
        a ^= (b ^ d ^ (a << 2) ^ (b << 18) ^ (c << 28)) & 0xffffffff
        state[index] = a
        self._index = index
        return a

    def random(self):
        #return self.getrandbits(53) * (1.0 / 9007199254740992.0)
        return self._genrand_int32() * 2.3283064365386963e-10

    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):
        self.__init__(n)

    def getrandbits(self, k=32):
        ret = self._genrand_int32()
        for _ in range((k - 1) // 32):
            ret = (ret << 32) ^ self._genrand_int32()
        return ret >> ((- k) % 32)

# EOF

メルセンヌ・ツイスタ（以下MTと略します）は、とても優秀なアルゴリズムです。
ただ、一般に出回っている初期化アルゴリズムは、seedが32ビットで設計されているため、その周期を生かし切れないケースが稀にあります。

例えばトランプのシャッフルです。52枚のカードをシャッフルするのに必要なパターンは、52の階乗（52×51×・・・×2×1)です。これは だいたい8.0×10の52乗になり、32ビットで表現できるパターン数をはるかに上回ります。

ちなみにPythonのrandomモジュールも初期化の種は32ビットの整数なので、同様の問題点があります。

乱数の周期は問題ないのですが、seedの桁数が不足しているんですよね。
かといってseedの桁数を増やしていくと、そのseedの生成にハッシュが必要になってしまいます。

・・・まてよハッシュ？
SHA512の出力を、そのまま内部状態の初期化に使えば良いんじゃないだろうか？

ただし、MTは、もっとも内部状態がすくないもので、521ビットの内部状態を持ちます。

・・・むむ、9ビット足りない。なんかすっきりしないなー。

そんなもやもやを抱えてWebを放浪していたところ、WELLにたどり着きました。

WELLとは？

• MTの開発者も関わっているMTより新しいアルゴリズム
• MTより少し計算コストがかかる
• MTより乱数の品質（高次の均等分布や零超過状態の回復）がより良い。
• サポートしている内部状態は512～444972ビット

・・・512ビット！ちょうどSHA512で初期化できるぞ！

ということで、WELL512＋SHA512でシャッフルすることにしました。

でも、バベルみたいにデッキが250枚くらいになると、512ビットでも不足するんですね。
そこまで考慮すると、WELL1024 + SHA512×2になりますね。でも、これはちょっと美しくない感じがするので、バベルは無視することにしました！まぁ人間のシャッフルよりはずっと良いですしね！