Xorshift の Python 実装

本日は

Wikipedia に記載されている Xorshift をPythonで実装しました：
//https://ja.wikipedia.org/wiki/Xorshift

Wikipediaの場合

//https://ja.wikipedia.org/wiki/Xorshift
# include <stdio.h>
# include <stdint.h>

uint32_t xor32(void)
{
    static uint32_t y = 2463534242;
    y = y ^ (y << 13);
    y = y ^ (y >> 17);
    y = y ^ (y << 5);
    return y;
}

uint32_t xor64(void)
{
    static uint64_t x = 88172645463325252ULL;
    x = x ^ (x << 13);
    x = x ^ (x >> 7);
    x = x ^ (x << 17);
    return x;
}

uint32_t xor96(void)
{
    static uint32_t x = 123456789;
    static uint32_t y = 362436069;
    static uint32_t z = 521288629;
    uint32_t t;
    t = (x ^ (x << 3)) ^ (y ^ (y >> 19)) ^ (z ^ (z << 6));
    x = y;
    y = z;
    z = t;
    return z;
}

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

int main(void)
{
    int count = 100;
    for (int i = 0; i < count; ++i)
    {
        printf("%u\n", xor32());
    }
}

clangなどでuint32_tなどを使うときは stdint.h をインクルードしないといけないようですね．

Python 版

Pythonだと整数値が任意の精度で計算されてしまうので

• XorShiftの実装 [C++/D言語/C#/VB/JS/Python/HSP/Lisp]
• 关于Xorshift随机数发生器算法

にならい0xFFFFFFFFでマスクをかけています．

　（Edit: numba jit が効果なかったので外しました．)

xorshift.py

"""
implement xorshift with python
https://ja.wikipedia.org/wiki/Xorshift
https://qiita.com/yosgspec/items/e4287262f8dbea2aa815
https://ask.helplib.com/1591921
"""

def xorshift(generator, seed=None):
    ret = seed
    def inner():
        nonlocal ret
        if ret is None:
            ret = generator()
        else:
            ret = generator(*ret)
        return ret[-1]
    return inner


def xor32(y=2463534242):
    y = y ^ (y << 13 & 0xFFFFFFFF)
    y = y ^ (y >> 17 & 0xFFFFFFFF)
    y = y ^ (y << 5 & 0xFFFFFFFF)
    return y & 0xFFFFFFFF,


# def xor64(_x=88172645463325252, x=88172645463325252):
#    _x = _x ^ (_x << 13)
#    _x = _x ^ (_x >> 7)
#    _x = _x ^ (_x << 17)
#    return _x, _x & 0xFFFFFFFF


def xor96(x=123456789, y=362436069, z=521288629):
    t = (x ^ (x << 3 & 0xFFFFFFFF)) ^ (
        y ^ (y >> 19 & 0xFFFFFFFF)) ^ (
        z ^ (z << 6 & 0xFFFFFFFF))
    x = y
    y = z
    z = t
    return x, y, z


def xor128(x=123456789, y=362436069, z=521288629, w=88675123):
    t = x ^ (x << 11) & 0xFFFFFFFF
    x = y
    y = z
    z = w
    w = (w ^ (w >> 19)) ^ (t ^ (t >> 8)) & 0xFFFFFFFF
    return x, y, z, w


def uniform(rand, begin=0, end=1):
    def inner():
        return begin+(rand())/(int(0xFFFFFFFF)/(end-begin))
    return inner


def calc_pi(generator):
    u01 = uniform(generator)
    counter = 0
    N = 100000
    for i in range(N):
        x = u01()
        y = u01()
        if x*x+y*y < 1.0:
            counter += 1
    print(4.0*counter/N)


def apply_example():
    random32 = xorshift(xor32)
    calc_pi(random32)
    #random64 = xorshift(xor64)
    #calc_pi(random64)
    random96 = xorshift(xor96)
    calc_pi(random96)
    random128 = xorshift(xor128)
    calc_pi(random128)

def main():
    random32=xorshift(xor32)
    for i in range(100):
        print(random32())

    apply_example()
if __name__ == '__main__':
    main()

xorshift 関数内にretの値を保存しておくことで呼び出しごとに異なる値を出すような仕組みを作っています． nonlocal の典型的な使い方ですね．

• あまり知られていないPythonの言語仕様（Python3.4以降対応）

できてないこと．

xorNN 関数でNNが３２，９６，１２８の場合はCのコードの結果と一致します．N=６４の場合はうまく一致できないのでコメントアウトしています．Cのほうで

static uint64_t x = 88172645463325252ULL;

の定義を変えることでPythonの実装と合わせることができますが，それはいかがなものかと・・・．