すごい乱数生成アルゴリズム「xorshift」

Python

みなさん、こんにちは、今回は乱数の話です。


特に複数機種でのコンシューマ機でゲームを開発をしていると、機種間で乱数値を統一するために乱数生成アルゴリズムを自作しますよね。


そこでよく使われるアルゴリズムが「線形合同法」です、内容は至って簡単で、以下の漸化式を使います。


X_{n+1} = \left( A \times X_n + B \right) \bmod M


A,B,Mは定数で、どの値が入るかは処理系依存です。
例えばUnixなどの処理系ではA=1103515245,B=12345,M=2147483647などが入ります。
C言語ですと以下のようになります。

static unsigned int x=1;

void srand(unsigned int s) { x=s; }

unsigned int rand() {
    x=x*1103515245UL+12345UL;
    return x&2147483647UL;
}


この「線形合同法」は計算が簡単で高速ですから、いろいろな環境と言語で使われています。
（C言語やJava、C#などの標準乱数は、このアルゴリズムが使われています）


しかし、線形合同法の生成する乱数は品質に非常に問題ありなんです。
周期はどう頑張っても2^32以上は期待できませんし、特に下位ビットの方は周期的な偏りがあり品質が非常に悪い。
最近では下位ビットを捨てて、品質を上げている処理系もあります。


品質を求めるなら「メルセンヌツイスター」でしょうね、(2^219937)-1とほぼ無限とも言える周期で、値もまんべんなく分布しています。
でも、計算時間が少し長いので、リアルタイムゲームではちょっと使いにくいのが玉にキズですね。
（でも最近は改良型のSFMTが出てきて、それも変わりつつあるんですが・・・）

簡単だけどすごいやつ「xorshift」

で、今回紹介するのが「xorshift」です、こいつがスゴイ！
メルセンヌツイスターには及ばないものの、(2^128)-1の周期を持っており値の偏りも少なく、品質が非常に高いのが特徴です。
しかも驚くことに、計算はビットシフトとXORだけでなので非常に高速です。
以下、Pythonでの実装例です。

#coding:utf8
x,y,z,w = (123456789,362436069,521288629,88675123)

def seed(s):
    global x,y,z,w
    x,y,z,w = s
    if x+y+z+w <= 0:
        raise ValueError, "Please do not substitute 0 for seed."

def randint(a,b):
    global x,y,z,w
    t = (x^(x<<11))
    x,y,z = y,z,w
    w = w=(w^(w>>19))^(t^(t>>8))
    return (a+w)%b

#乱数のシードは4つの整数値を指定
seed([12345678,4185243,776511,45411])

#乱数生成
for _ in range(100):
    print randint(0,100)


こんな簡単なコードなのにちゃんと乱数してますね。
計算もXORとシフトだけなので、計算速度は線形合同法と比べても遜色はありません、品質は断然xorshiftが上ですけどね。
xorshiftの論文は、以下のページからPDFファイルをダウンロードできます。（英文）


http://www.jstatsoft.org/v08/i14/


論文には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の標準乱数は、C言語で書かれたメルセンヌツイスターです。なので、Pythonでxorshiftを実装する意味はあんまないです（ならなぜ紹介したっ！）

何が言いたいかといいますと、Python最高！！