乱数生成アルゴリズム Xorshift を random.Random のサブクラスとして実装してみた

すごい乱数生成アルゴリズム「xorshift」 - Pashango’s Blog の Python 最高！！につられて。 random.Random を継承する形で Xorshift アルゴリズムを。なぜつくったのかまったくの謎だけど Python 最高！！だからしかたない。

メルセンヌ・ツイスタ (Mersenne Twister) のほうが高品質なので、出番はたぶんない。ホントなぜつくった、オレ？

random.Random のサブクラスにすると randint や choice などのメソッドを動かすことができるのがメリット。

Python 版

まずは Pure Python で。

核となる乱数生成部は _genrand_int32 メソッド。名前の通り xor とシフト演算しかないのと、状況保持が unsinged long 4 つだけでできるのが特徴。論文にあったオリジナルの C コードは以下。

unsigned long xor128(){ 
    static unsigned long x=123456789,y=362436069,z=521288629,w=88675123; 
    unsigned long t; 
    t=(x^(x<<11));x=y;y=z;z=w; return( w=(w^(w>>19))^(t^(t>>8)) ); 
}

これを Python 用に書き換えたものが _genrand_int32。値が 32 ビット以上にならないように、左シフトしたところは & 0xffffffff している。ちなみに、オリジナルのコードにでてくる 4 つの値 x, y, z, w は、属性 _state にリストで保持している。

seed メソッドは http://unkar.jp/read/pc11.2ch.net/tech/1192628099#l76 を、その他は Python 2.6.2 の Lib/random.py と Modules/_randommodule.c を参考に作成。参考ではなくパクりかも。

なお、 jumpahead と getrandbits メソッドは未実装。

xorshift_py.py

# coding: utf-8
import random
from os import urandom as _urandom
from binascii import hexlify as _hexlify

class Xorshift(random.Random):

    VERSION = (u'Xorshift', 1) # used by getstate/setstate

    def __init__(self, x=None):
        """Initialize an instance.

        Optional argument x controls seeding, as for seed().
        """
        self._state = [1812433254, 3713160357, 3109174145, 64984499]
        self.seed(x)
        self.gauss_next = None

    def _genrand_int32(self):
        t = self._state[0] ^ ((self._state[0] << 11) & 0xffffffff)
        self._state[:] = self._state[1], self._state[2], self._state[3], \
                (self._state[3] ^ (self._state[3] >> 19)) ^ (t ^ (t >> 8))
        return self._state[3]

    def random(self):
        """Get the next random number in the range [0.0, 1.0)."""
        a = self._genrand_int32() >> 5
        b = self._genrand_int32() >> 6
        return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0)

    def seed(self, a=None):
        """Initialize internal state from hashable object.

        None or no argument seeds from current time or from an operating
        system specific randomness source if available.

        If a is not None or an int or long, hash(a) is used instead.
        """
        if a is None:
            try:
                a = long(_hexlify(_urandom(16)), 16)
            except NotImplementedError:
                import time
                a = long(time.time() * 256) # use fractional seconds
        elif isinstance(a, (int, long)):
            a = abs(a)
        else:
            a = hash(a) & 0xffffffff

        for i in range(1, 5):
            a = (1812433253 * (a ^ (a >> 30)) + i) & 0xffffffff
            self._state[i-1] = a

        self.gauss_next = None

    def getstate(self):
        """Return internal state; can be passed to setstate() later."""
        return self.VERSION, tuple(self._state), self.gauss_next

    def setstate(self, state):
        """Restore internal state from object returned by getstate()."""
        version = state[0]
        if version == (u'Xorshift', 1):
            self._state = list(state[1])
            self.gauss_next = state[2]
        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):
        raise NotImplementedError

def main():
    r = Xorshift(1000)
    for _ in range(5):
        print '%.6f' % r.random()
    for _ in range(5):
        print '%08d' % r.randint(0, 100000000)
    s = 'abcdefg'
    for _ in range(5):
        print '%s' % r.choice(s)

if __name__ == '__main__':
    main()

C 版

C 言語で Python Extension を作成。 Python/C API の練習がしてみたかった、ただそれだけ。乱数生成器なんだから速度は重要だろう、などといった高尚な動機はなかったりする。 Pure Python 版と同じ動作をする、はず。

Python で xorshift.Xorshift クラスをつくっておいて、状態保持と乱数生成を C で書いた _xorshift.Xorshift クラスに委譲。

最初は多重継承で作ってみていたのだけれども C の自作クラスと random.Random を多重継承したら TypeError だった。 Python/C API 暦 10 時間未満の自分では対処できず。

TypeError: Error when calling the metaclass bases
    multiple bases have instance lay-out conflict

乱数生成に用いる 4 つの値 x, y, z, w は XorshiftObject 構造体の unsigned long の state[4] 配列に保持。 Python オブジェクトには変換せず、 C の世界に置きっぱなしにすることで速度アップを図る。

しかし Python 側からもこれらの値を操作できるようにしておかないと getstate, setstate メソッドの作成に支障をきたすので、これらの入出力を行う xorshift_getstate, xorshift_setstate 関数を作成し、メソッドとして公開しておく。

それでも自分の環境では random.Random より遅かった。コンパイラオプションとか、コードのこまかい調整で早くなるのかもしれないけれど。

環境

Windows XP Professional
Microsoft Visual C++ 2008 Express Edition
Python 2.6.2

試していたのは Windows でだったけれども、特に変わったことはしていないので Python が動く環境なら OS 問わず setup.py build できる気がする。

むしろ Windows のほうが C Extension を作るの、難しいような。 Python 2.5 の頃は Visual Studio 2003 を要求されて、用意できなくて mingw32 使ったりした思い出が。とはいえ Python 2.6 は 2008 でいいからだいぶ敷居がさがっている、はず。

setup.py

from distutils.core import setup, Extension

py_modules = ['xorshift']
ext_modules = [Extension('_xorshift', sources = ['_xorshift.c'])]

setup (name = 'xorshift',
       version = '0.1',
       author='fgshun',
       author_email='fgshun@lazy-moon.jp',
       url='http://d.hatena.ne.jp/fgshun/',
       py_modules=py_modules,
       ext_modules=ext_modules)

_xorshift.c

#include "Python.h"
#include <time.h>

typedef struct {
    PyObject_HEAD
    unsigned long state[4];
} XorshiftObject;

static PyTypeObject Xorshift_Type;

static unsigned long
xorshift_genrand_int32(XorshiftObject *self)
{
    unsigned long t; 

    t = self->state[0] ^ (self->state[0] << 11);
    self->state[0] = self->state[1];
    self->state[1] = self->state[2];
    self->state[2] = self->state[3];
    self->state[3] = (self->state[3] ^ (self->state[3] >> 19)) ^
                     (t ^ (t >> 8)); 
    return self->state[3];
} 

static PyObject *
xorshift_random(XorshiftObject *self)
{
    unsigned long a, b;
    a = xorshift_genrand_int32(self) >> 5,
    b = xorshift_genrand_int32(self) >> 6;
    return PyFloat_FromDouble(
            (a * 67108864.0 + b) * (1.0 / 9007199254740992.0));
}

static PyObject *
xorshift_seed(XorshiftObject *self, PyObject *args)
{
    PyObject *arg = NULL;
    PyObject *n = NULL;
    unsigned long seed;
    int i;
    unsigned long *state;

    if (! PyArg_UnpackTuple(args, "seed", 0, 1, &arg))
        return NULL;

    if (arg == NULL || arg == Py_None) {
        time_t now;
        time(&now);
        seed = (unsigned long)now;
    }
    else if (PyInt_Check(arg) || PyLong_Check(arg)) {
        n = PyNumber_Absolute(arg);
        if (n == NULL)
            return NULL;
        seed = PyInt_AsUnsignedLongMask(n);
        Py_DECREF(n);
    }
    else {
        long hash = PyObject_Hash(arg);
        if (hash == -1)
            return NULL;
        seed = (unsigned long)hash;
    }
    
    state = self->state;
    seed = seed & 0xffffffffUL;
    for (i=1; i<4+1; i++) {
        state[i-1] = seed = 
            (1812433253UL * (seed ^ seed >> 30)) + i;
        seed &= 0xffffffffUL;
        state[i-1] = seed;
    }

    Py_INCREF(Py_None);
    return Py_None;
}

static PyObject *
xorshift_getstate(XorshiftObject *self)
{
    PyObject *state = NULL;

    state = Py_BuildValue("kkkk",
            self->state[0], self->state[1], self->state[2], self->state[3]);
    if (state == NULL)
        return NULL;
    return state;
}

static PyObject *
xorshift_setstate(XorshiftObject *self, PyObject *args)
{
    unsigned long x, y, z, w;

    if (! PyArg_ParseTuple(args, "(kkkk)", &x, &y, &z, &w))
        return NULL;
    self->state[0] = x;
    self->state[1] = y;
    self->state[2] = z;
    self->state[3] = w;

    Py_INCREF(Py_None);
    return Py_None;
}

static PyObject *
xorshift_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
{
    XorshiftObject *self;
    PyObject *tmp;

    self = (XorshiftObject *)type->tp_alloc(type, 0);
    if (self == NULL) 
        return NULL;
    tmp = xorshift_seed(self, args);
    if (tmp == NULL) {
        Py_DECREF(self);
        return NULL;
    }
    Py_DECREF(tmp);
    return (PyObject *)self;
}

static PyMethodDef xorshift_methods[] = {
    {"random", (PyCFunction)xorshift_random, METH_NOARGS,
        PyDoc_STR("random() -> x in the interval [0, 1).")},
    {"seed", (PyCFunction)xorshift_seed, METH_VARARGS,
        PyDoc_STR("seed([n]) -> None.")},
    {"getstate", (PyCFunction)xorshift_getstate, METH_NOARGS,
        PyDoc_STR("getstate() -> tuple containging the current state.")},
    {"setstate", (PyCFunction)xorshift_setstate, METH_VARARGS,
        PyDoc_STR("setstate() -> None. Restores generator state.")},
    {NULL} /* sentinel */
};

PyDoc_STRVAR(xorshift_doc, 
        "create a random number generator with its own internal state.");

static PyTypeObject Xorshift_Type = {
    PyVarObject_HEAD_INIT(NULL, 0)
    "_xorshift.Xorshift",        /*tp_name*/
    sizeof(XorshiftObject),      /*tp_basicsize*/
    0,                           /*tp_itemsize*/
    /* methods */
    0,                /*tp_dealloc*/
    0,                /*tp_print*/
    0,                /*tp_getattr*/
    0,                /*tp_setattr*/
    0,                /*tp_compare*/
    0,                /*tp_repr*/
    0,                /*tp_as_number*/
    0,                /*tp_as_sequence*/
    0,                /*tp_as_mapping*/
    0,                /*tp_hash*/
    0,                /*tp_call*/
    0,                /*tp_str*/
    PyObject_GenericGetAttr,    /*tp_getattro*/
    0,                /*tp_setattro*/
    0,                /*tp_as_buffer*/
    Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE,    /*tp_flags*/
    xorshift_doc,     /*tp_doc*/
    0,                /*tp_traverse*/
    0,                /*tp_clear*/
    0,                /*tp_richcompare*/
    0,                /*tp_weaklistoffset*/
    0,                /*tp_iter*/
    0,                /*tp_iternext*/
    xorshift_methods, /*tp_methods*/
    0,                /*tp_members*/
    0,                /*tp_getset*/
    0,                /*tp_base*/
    0,                /*tp_dict*/
    0,                /*tp_descr_get*/
    0,                /*tp_descr_set*/
    0,                /*tp_dictoffset*/
    0,                /*tp_init*/
    0,                /*tp_alloc*/
    xorshift_new,     /*tp_new*/
    _PyObject_Del,    /*tp_free*/
    0,                /*tp_is_gc*/
};

static PyMethodDef module_methods[] = {
    {NULL} /* Sentinel */
};

PyDoc_STRVAR(module_doc,
"The Xorshift random number generator.");

PyMODINIT_FUNC
init_xorshift(void)
{
    PyObject *m;

    if (PyType_Ready(&Xorshift_Type) < 0)
        return;

    m = Py_InitModule3("_xorshift", module_methods,
            module_doc);

    if (m == NULL)
        return;

    Py_INCREF(&Xorshift_Type);
    PyModule_AddObject(m, "Xorshift", (PyObject *)&Xorshift_Type);
}

xorshift.py

# coding: utf-8
import random
from os import urandom as _urandom
from binascii import hexlify as _hexlify
import _xorshift

class Xorshift(random.Random):

    VERSION = (u'Xorshift', 1) # used by getstate/setstate

    def __init__(self, x=None):
        """Initialize an instance.

        Optional argument x controls seeding, as for seed().
        """
        self._xorshift = _xorshift.Xorshift(x)
        self.gauss_next = None

    def random(self):
        """Get the next random number in the range [0.0, 1.0)."""
        return self._xorshift.random()

    def seed(self, a=None):
        """Initialize internal state from hashable object.

        None or no argument seeds from current time or from an operating
        system specific randomness source if available.

        If a is not None or an int or long, hash(a) is used instead.
        """
        if a is None:
            try:
                a = long(_hexlify(_urandom(16)), 16)
            except NotImplementedError:
                import time
                a = long(time.time() * 256) # use fractional seconds

        self._xorshift.seed(a)
        self.gauss_next = None

    def getstate(self):
        """Return internal state; can be passed to setstate() later."""
        return self.VERSION, self._xorshift.getstate(), self.gauss_next

    def setstate(self, state):
        """Restore internal state from object returned by getstate()."""
        version = state[0]
        if version == (u'Xorshift', 1):
            self._xorshift.setstate(state[1])
            self.gauss_next = state[2]
        else:
            raise ValueError("state with version %s passed to "
                             "setstate() of version %s" %
                             (version, self.VERSION[0]))

    def jumpahead(self, n):
        raise NotImplementedError

    def getrandbits(self, k):
        raise NotImplementedError

def main():
    r = Xorshift(1000)
    for _ in range(5):
        print '%.6f' % r.random()
    for _ in range(5):
        print '%08d' % r.randint(0, 100000000)
    s = 'abcdefg'
    for _ in range(5):
        print '%s' % r.choice(s)

if __name__ == '__main__':
    main()

銀月の符号

Python 使い見習いの日記・雑記