トップ
乱数ライブラリ

良い乱数・悪い乱数

C言語標準ライブラリの乱数rand( )は質に問題があり、禁止している学会もある。 他にも乱数には様々なアルゴリズムがあるが、多くのものが問題を持っている。 最も多くの人に使われている乱数であろう Visual Basic の Rnd の質は最低である。

そもそも乱数とは

乱数とは、本来サイコロを振って出る目から得られるような数を意味する。 このような乱数は予測不能なものである。 しかし、計算機を使って乱数を発生させた場合、 次に出る数は完全に決まっているので、予測不能とはいえない。 そこで、計算機で作り出される乱数を疑似乱数(PRNG)と呼び区別することがある。 ここでは、特にことわらない限り乱数とは疑似乱数のことを指すとする。

計算機でソフト的に乱数を発生させることの最大のメリットは、 再現性があることである。 初期状態が同じであれば、発生する乱数も全く同じものが得られる。 このことは、パラメータを変更して実験する場合、 乱数による揺らぎを完全に無くすのに便利である。 もちろん、サイコロを振って得られる乱数には再現性がない。

ほとんどの乱数生成プログラムは、種と呼ばれる数値を受け取って、 それに対応する乱数列を発生する。 実験のとき、種を別のものにすることによって、 簡単に異なる乱数列を得ることができる。

「乱数はもともとでたらめなので、いい加減な方法を使っても問題は起きない」 という考えがある。この考えは、危険である。 例えば、現在C言語の入門書などで広く紹介されている rand( ) という乱数は、 期待するほどでたらめではなく、ある意味で規則的な数列を返す。 従って、例えば rand( ) を使って実験を行ない、 統計処理した場合、結果の信頼性が低くなる。

このページの目的は、従来の乱数の問題点を指摘し、 すばらしい乱数、 「メルセンヌツイスタ」 の普及を促進することである。

また、乱数のことについて深く解説した杉田洋氏のページ モンテカルロ法,乱数,および疑似乱数も紹介しておく。

従来の乱数生成プログラムの問題点

rand( )その1

UNIX の /usr/ucb/cc を解析した結果、 rand( ) は以下の動作をすることがわかった。

static long x=1;
void srand(long s) { x=s; }
long rand() { x=x*1103515245+12345; return x&2147483647; }

これは、非常にシンプルな線形合同法である。 この乱数の最下位ビットは0と1の繰り返しになる。 すなわち、偶数と奇数が交互に生成される。 このことから、この乱数で下位ビットを乱数として使うのは危険であることがわかる。 また、ある乱数が得られたら、次に現れる乱数が1種類しかないという欠点も持つ。

rand( )その2

多くの rand( ) が以下のアルゴリズムによって生成されている。

static long x=S;
int rand() { x=x*A+C; return(int)(x>>16)&32767; }
void srand(long s) { x=s; if (F) rand(); }

/usr/5bin/cc, LSI-C  A=1103515245, C=12345,   F=0, S=1
TurboC 1.5           A=22695477,   C=1,       F=0, S=1
Visual C++           A=214013,     C=2531011, F=0, S=1。
Borland C++          A=22695477,   C=1,       F=1, S=22695478。

これらのソースも筆者が独自に解析した結果、得られたものである。

rand( )その1、その2に共通していることは、線形合同法を使っていることである。 すなわち、変数xは漸化式

X[n]=A*X[n-1]+C (mod 4294967296)

に基づいて計算される。 (mod 4294967296)とは4294967296で割った余りという意味である。 多くの計算機では余りを求める演算をわざわざ行なわなくても、 掛け算と足し算の結果オーバーフローした部分を無視することによって、 自動的にこの結果が得られる。

X[n]の持つ重要な性質は、最下位ビットの周期が2、すなわち0と1が交互に現れ、 その隣のビットの周期が4、さらにその隣のビットは8、となり、 最上位ビットの周期は4294967296となることである。 これらのrand( )はxの最上位ビットを無視するので、全体の周期は2147483648になる。 rand( )のマニュアルには周期は4294967296と書いてあるがそれは間違いである。 実際、この周期は、乱数を使った実験では短すぎる場合が多い。 また、以前パチスロで一定周期で当たりが出やすくなるという問題が起きたが、 その原因は rand( ) を使っていたためと私は推測する。

rand( )その2 はX[n]の下位16ビットを捨て、 上位16ビットを右シフトしてそのうちの下位15ビットを返す。 従って最下位ビットでも周期は131072になるので、 「rand( ) その1」のような下位ビットを使った極端な問題は起きない。 しかし、乱数の返す数値の範囲が、0から32767までしかない。 言い換えれば、精度が15ビットしかない。 そのため、例えば、Box-Muller法で正規分布に従う乱数を発生させたり、 逆関数法で指数分布に従う乱数を発生させたりするとき、 端の方で理論的な分布と、生成された乱数の分布との違いが大きくなる。 下の表は全周期で最大値や最小値を調べた結果である。

              正規分布            指数分布
              最小値   最大値     最大値
rand()その1  -6.53σ  6.35σ     21.49λ
rand()その2  -4.56σ  4.56σ     10.40λ

rand( )の一番大きな欠点は、移植性がないことである。 rand( )その1の最大値は 2147483647 であるのに対し、 rand( ) その2は 32767 なので、異なるコンパイラでプログラムの移植するときに、 ソースコードを書き換えなくてはならない。 また、rand( )その2の間でもパラメータが違うものがあるので、 一方のコンパイラで起きたことを、他方のコンパイラで再現することが不可能である。

さらに、rand( )に共通する欠点として、違う種をsrand( )に与えても、 同じ乱数列が発生する場合があることも挙げられる。 例えば、srand(0)とsrand(0x80000000)で初期化された乱数列は、互いに同じものになる。

Visual Basic の Rnd

Visual Studio やエクセルなどで使われる Visual Basic の Rnd は簡略化して書けば以下のようになっている。

static long x=327680;

float Rnd(void)
{
    x=x*16598013+12820163&16777215;
    return x*(1.0/16777216.0);
}

見てわかるように、24ビット線形合同法を使っている。 周期はたったの 1677万しかない。ちょっとした実験なら一周してしまううえ、精度も24ビットしかない。

もっと大きな問題は、乱数を初期化する Randomize にある。 これは与えられた種を16ビットの整数に変換して初期化する。従って、65536種類の 系列しか得られない。違う種を指定して実験しても、 かなりの確率で、まったく同じ系列が返される危険がある。

参考までに、Rnd と Randomize と全く同じ動作をする Rnd2 と Randomize2 を以下に示す。

自前のRnd ( http://www001.upp.so-net.ne.jp/isaku/rnd.html )

C言語やVisual Studioなら乱数のプログラムを書いて、対処できるが、 エクセルの場合、整数のオーバーフローを認めていないので、性能のよい乱数を使いたいなら、 C言語などでダイナミック・リンク・ライブラリのプログラムを書いて、それを呼び出すしかない。 多くの乱数は 32 ビット整数演算をしているのでエクセルVBAの Long 型では計算できない。 以下のプログラムは後述する Xorshift とメルセンヌツイスタを 31 ビットにしてVBAで書いたもである。

VBエディターを開いて、標準モジュールを追加し、そこに貼り付ければ使うことができる。 Xorshift は周期を最大化するパラメータを探すプログラムを自分で書いて 31×4 ビットのものを得た。 メルセンヌツイスタは、 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html にある、Dynamic Creatorを使って 31 ビットに対応したパラメータを計算して求めた。 メルセンヌツイスタの方が品質も良く周期も長いが、Xorshift に比べて遅い。Core 2 Duo E6600の エクセル2007で1億個発生させるのに Xorshift で約1分、メルセンヌツイスタで約5分かかった。 スピードを求めるなら乱数ライブラリのダイナミック・リンク・ライブラリを使うべきであろう。 (メルセンヌツイスタの周期:4.31e6001、Xorshiftの周期:2.12e32)

drand48( )

UNIXでは、drand48( )と呼ばれる乱数パッケージを使うことができる。 使い方は、man drand48 で知ることができる。 以下にdrand48( )の一部の機能を実現するソースを示す。

static long long  x=0x1234ABCD330E;

double drand48() /* 0.0以上1.0未満 */
{
    x = x*0x5DEECE66D + 0xB;
    return (x&0xFFFFFFFFFFFF)*(1.0/281474976710656.0);
}

long lrand48() /* 0以上2147483647以下 */
{
    x = x*0x5DEECE66D + 0xB;
    return(long)(x>>17)&0x7FFFFFFF;
}

long mrand48() /* -2147483648以上214748367以下 */
{
    x = x*0x5DEECE66D + 0xB;
    return(long)(x>>16)&0xFFFFFFFF;
}

void srand48(long s)
{
    x = s;
    x = (x<<16)+0x330E;
}

ソースを見てわかるように、drand48( )も線形合同法を使っている。 rand( )は32ビットの線形演算を行い、そのうち31ビットを使っていた。 これに対して、drand48( )は64ビットの線形演算を行い、そのうち、48ビットを使っている。 従ってdrand48( )は2の48乗の周期を持つ。 この周期の長さは rand( ) に比べて十分長い。 ただし、mrand48( )の最下位ビットの周期は 131072 しかないことに注意すること。 また、64ビットの掛け算を行っているため、 実行時間が rand( ) に比べかなり遅くなっている。

drand48( )の長所として、移植性の高さが挙げられる。 ほとんどの UNIX で drand48( ) が使用でき、全く同じ動作をする。 Visual C++で drand48( ) が使用できないのが残念に思われる。

線形合同法の欠点

rand( ) その1、rand( ) その2、drand48( ) は線形合同法で乱数を 発生させている。この方法に共通する欠点は、ある乱数が得られたら、 次に現れる乱数が限られてしまうことだ。とくに、rand( ) その1と、 drand48( ) は、次の乱数が完全に決まってしまう。 この欠点は乱数を組にして使うと問題になる。 例えば2次元ベクトル(x,y)を線形合同法で発生させると、 限られたベクトルしか現れない。

ちなみに、メルセンヌツイスタでは623次元のベクトルを発生させても、 全周期で、全ての組み合わせのベクトルが現れる。

random( )

多くの UNIX では random( ) という乱数が使える。 使い方は man 3 random で知ることができる。 以下に random( ) を簡略化したソースを示す。
unsigned x[32];

int random(void)
{
    unsigned n=(unsigned)x[0];

    if (n==31) n=1; else n++;
    x[0]=n; x[n]+=(n>3?x[n-3]:x[n+31-3]);
    return(x[n]>>1)&0x7FFFFFFF;
}

/* 以前の srandom */
void srandom(int s)
{
    int i;

    x[0]=3; x[1]=(unsigned)s;
    for(i=2;i<=31;i++) x[i] = 1103515245*x[i-1] + 12345;
    for(i=0;i<310;i++) random();
}

/* 現在の srandom */
void srandom(int seed) {
    int i; long long s=seed;
    if (s==0) s=1;
    x[0]=3; x[1]=(unsigned)s;
    for (i=2;i<=31;i++) {
        s=(16807*s)%2147483647;
        if (s<0) s+=2147483647;
        x[i]=(unsigned)s;
    }
    for (i=0;i<310;i++) random();
}

random( )は足し算法を使っている。すなわち、変数Xnは漸化式

X[n]=X[n-3]+X[n-31] (mod 4294967296)

に基づいて計算される。最下位ビットの周期は2147483648で、 全体の周期は約2の62乗である。 random( )の優れた点は生成するときの実行時間が短いことである。 乱数の質が問題にならず、速度だけが問題になる時は、この乱数を使うとよい。

しかし、random( ) には統計的な欠陥が見つかっているので注意が必要である。 中心極限定理により、乱数を多数加えると正規分布に従うはずであるが、 random( ) を使った場合は正規分布から少し外れてしまう。

また、新しい srandom() の場合、種を0にした場合と1にした場合に、 全く同じ結果となる。この事は man 3 random、man 3 random_r と得られるマニュアルに 書かれておらず、ソースコードのコメントにも書かれていない。マルチスレッドで、 srandom_r( ) にスレッドの番号を渡すようなプログラムを組むような場合注意が必要である。

最下位ビットの問題

random( )のマニュアルには全てのビットがランダムで、 使用可能と書かれているが、以前の srandom( ) の初期化がまずいため、 下位ビットを使うのが好ましくないことがある。

例えば、以下のプログラムは、以前の random( ) を使って種を4つおきにして初期化し、 62 個の最下位ビットを表示するものだが、

/* TEST1 */
main(void)
 {
    int i,j,k;
    for (i=0;i<4;i++) for (j=0;j<4;j++) {
       printf("%4d ",i+j*4 );
       srandom(i+j*4);
       for (k=0;k<62;k++) printf("%ld",random()&1);
       printf("\n");
   }
   return 0;
}

結果は、

   0 11010110001011101100111001101000101011011001111010001101111001
   4 11010110001011101100111001101000101011011001111010001101111001
   8 11010110001011101100111001101000101011011001111010001101111001
  12 11010110001011101100111001101000101011011001111010001101111001
   1 10101101101110111000010011110101100110000111101101100001101111
   5 10101101101110111000010011110101100110000111101101100001101111
   9 10101101101110111000010011110101100110000111101101100001101111
  13 10101101101110111000010011110101100110000111101101100001101111
   2 01110010110101000011011001011011100000011110010000101001000111
   6 01110010110101000011011001011011100000011110010000101001000111
  10 01110010110101000011011001011011100000011110010000101001000111
  14 01110010110101000011011001011011100000011110010000101001000111
   3 00001001010000010111110011000110101101000000000111000101010001
   7 00001001010000010111110011000110101101000000000111000101010001
  11 00001001010000010111110011000110101101000000000111000101010001
  15 00001001010000010111110011000110101101000000000111000101010001

となる。以前の random( ) では4種類の全く同じ系列が順番に現れることがわかる。 この問題は現在の random( ) では解決されている。 しかし、この問題はlrand48( )で生じる。 従って、種を変えて100回実験をしても、得られる結果は4通りしかない。 次の表に、各乱数の最下位ビットの特性をまとめてみた。

                周期       系列の種類の数
rand()その1    2          2
rand()その2    131072     131072
random()        4294967294 4
lrand48()       262144     4
mrand48()       131072     2

これを見てわかることは、最下位ビットを使って実験するとき、 以前の random( ) や lrand48( ) を使うよりも rand( )その2の方が、 ましな結果が得られるということである。

.NET の Random クラス

Random クラスは Knuth の 引き算法を使っている。基本的には random( ) と変わりがない。 以下にコンストラクタと Next( ) と NextDouble( ) を簡略化したソースを示す。 Next( ) の最大値が 2147483646 であることに注意

static int index,x[56];

void Init(int seed) {
   int i,j,k;

   seed=161803398-(seed<0?-seed:seed); x[55]=seed; k=1;
   for (i=1;i<55;i++) {
       j=(21*i)%55; x[j]=k; k=seed-k;
       if (k<0) k+=2147483647;
       seed=x[j];
   }
   for (j=1;j<5;j++) for (i=1;i<56;i++) {
       x[i]-=x[1+(i+30)%55];
       if (x[i]<0) x[i]+=2147483647;
   }
   index=0;
}

int Next(void) {
    int r;

    if (++index>=56) index=1;
    r=x[index]-x[index<=34?index+21:index-34];
    if (r<0) r+=2147483647;
    x[index]=r; return r;
}

double NextDouble(void) { return Next()*(1.0/2147483647.0); }

上位ビットの問題

例えば、以下のC#プログラムは、0 から 63 までを種として初期化し、 100 個の乱数の最上位ビットを表示するものだが、

/* TEST2 */
using System;
class Program {
    static int SkippedRand(int seed,int NumOfSkip) {
        Random r=new Random(seed);
        for (int i=0;i<NumOfSkip;i++) r.Next();
        return r.Next();
    }
    static void Main(string[] args) {
        for (int i=1;i<=100;i++) {
            Console.Write("{0,2} 番目:",i);
            for (int s=0;s<64;s++)
                Console.Write("{0}",SkippedRand(s,i-1)>>30);
            Console.WriteLine();
        }
    }
}

結果は、

 1 番目:1010101010101101010101010101010101001010101010101010101011010101
 2 番目:1001101100100110110011011001001101100110110010011001001101100100
 3 番目:1001101100100110110010011011001001101100100110110010011011001001
 4 番目:1110011000110011100110001100111001100011001110011000110011100110
 5 番目:0101001010101011010101010010101010100101010101101010101001010101
 6 番目:1000011110000111100001111000011110000111100001111000011100001111
 7 番目:1010101001010101011010101010010101010010101010110101010110101010
 8 番目:0101010101010101010101010101010101001010101010101010101010101010
 9 番目:1000011110000111100001111000011110000111100001111000011110000111
10 番目:0100100101101001011010010010110100101101001001011010010110100101
11 番目:0011011001100110011001001100110011001100100110011001100110010011
12 番目:0001100110001100111001100111001100011001100011001110011001110011
13 番目:1001001001101101100100100110110110010010011011011001001001101101
14 番目:0110101010101010101010101011010101010101010101010100101010101010
15 番目:1100100110110010011011001001101100100110110010011011001001101100
16 番目:0101101001011010010110100101101001011010010110100101101001011010
17 番目:1010110101001010100101011010100101011010101101010010101101010110
18 番目:1101001011010010110100101001011010010110100101101011010010110100
19 番目:1111111111110000000000000000000111111111111111111000000000000000
20 番目:0100101001011010110100101101011010010100101101011010010110101101
〜中略〜
65 番目:0000000000000000000000000000000000000000000000000000000000000000
〜以下省略〜

となり、長さがまちまちであるが、周期的になっている。 この性質が引起こす問題は、例えば種が小さいときは、 65番目の乱数は必ず小さい数になるといった規則性を持つことである。

この現象はrand( )、lrand48( )、mrand48( )、以前の random( ) でも起こる。 線形合同法と引き算法は根本的に異なるアルゴリズムであるが、 この現象は共通して起こる。 例えば、Visual C++のrand( )の場合最初の乱数の最上位ビットは、 種が-11〜5005の間で0である。また、Borland C++の場合、79番目の 乱数は7407〜7408の間隔で最上位ビットが不変である

さらに、この現象は種を一定間隔に離しても起こる。 以下のプログラムは種を100おきに変化させながら最上位ビットを表示する。

/* TEST3 */
using System;
class Program {
    static int SkippedRand(int seed,int NumOfSkip) {
        Random r=new Random(seed);
        for (int i=0;i<NumOfSkip;i++) r.Next();
        return r.Next();
    }
    static void Main(string[] args) {
        for (int i=1;i<=20;i++) {
            Console.Write("{0,2} 番目:",i);
            for (int s=0;s<64;s++)
                Console.Write("{0}",SkippedRand(s*100,i-1)>>30);
            Console.WriteLine();
        }
    }
}

 1 番目:1100110011001100110011001100110001100110011001100110011001100110
 2 番目:1011011011011011011011010010010010010010010110110110110110110100
 3 番目:1110000011111000001111100000111110000011111000001111100001111100
 4 番目:1101101101001001001001001011011011011011010010010010010010110110
 5 番目:0011110001111000111000011100011110001110000111000011100011110001
 6 番目:1101011010010100101101011010010100101101011010010100101101011010
 7 番目:1110011100111001110001100011000110001110011100111001110001100011
 8 番目:0111000111000111000111000111000110001110001110001110001110001110
 9 番目:1010010110100101101001011010010010110100101101001011010010110100
10 番目:0001111000011110000111100001111000011110000111100001111000011110
11 番目:0110110010010010010010010010010011011011011011011011011011001001
12 番目:0111111100000001111111000000011111110000000111111100000001111111
13 番目:1001100110011001100110011001001100110011001100110011001100110011
14 番目:0011111111111111111110000000000000000000011111111111111111110000
15 番目:1111111111111111111111111111111111111111111111111111111111111111
16 番目:0010101010100101010101101010101001010101010010101010110101010101
17 番目:1111000000000000001111111111111100000000000000111111111111110000
18 番目:1111111111111111111111111111111111111111111111111111111111111111
19 番目:1001101100100110110011011001001100100110110010011001001101100110
20 番目:0010011001101100110010011001001100110110011001001100110110011011

今度は15番目と18番目の乱数で最上位ビットの周期が長くなっている。

この現象の原因は解明したが、その説明はここでは省く。 ここ(英文) に詳しい解説がある。 この現象のため、以前の random( )、rand( )、lrand48( )、mrand48( )、drand48( )、 .NET の Random クラスを使って実験をする場合、種は慎重に選ばなければならない。 連続して実験を行う場合に、安易に、システム時間等を使って種を設定すると、 種が一定間隔になり、実験結果に悪影響を及ぼす可能性がある。種を選ぶときには、 例えば 123456789、987654321、12121212などのように、十分大きな桁で、 無関係になるように工夫しなければならない。

悪い乱数のリスクを視覚的に明らかにするも参照のこと。

一方、現在の random( ) やメルセンヌツイスタでは、このような心配はいらない。 違う種を選べば、相関のない系列が得られる。

JAVAの乱数

JAVA の乱数は48ビット線形合同法を使っている。そして、統計的に問題がある。 nextDouble() で diehard test(http://i.cs.hku.hk/~diehard/) を行ったところ、OQSO test と DNA test で異常な結果が起こったからだ。

さらに、初期化ミスのため、例えば以下のプログラムを JAVA で実行すると、

    public class RandTest
    {
        static java.util.Random r=new java.util.Random();
        public static void main (String[] args)
        {
            double x,Min=1,Max=0;
            for (int i=0;i<1024;i++) {
                r.setSeed(i); x=r.nextDouble();
                if (x>Max) Max=x;
                if (x<Min) Min=x;
            }
            System.out.println("Min="+Min+" Max="+Max);
        }
    }

結果は、

    Min=0.6753377750582709 Max=0.7669794809891215

となり、0から1023までの種で初期化したにもかかわらず、 最初の乱数は0.7近辺しか現われない。

また、Math.random( ) は種を与える機能がない。 従って、実験などを再現することができないので、使い物にならない。

そもそも JAVA は遅く、乱数を使ったシミュレーションに向いていない。 さらに、48ビット線形合同法を使っていることが、その遅さにわをかけている。

メルセンヌツイスタ

繰り返すがメルセンヌツイスタ http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html は、乱数としてすぐれている。

メルセンヌツイスタはM系列乱数(GFSR)の子孫である。M系列乱数とは、 組み合わせ論から導き出された乱数で、各ビットが独立してk次均等分布していた。 k次均等分布とはk次元のベクトルを取ったときにk個がすべて0の場合を除いて、 全ての組み合わせが全周期で等しい回数だけ現れるというものである。

GFSRは各ビットが独立して生成されていた。メルセンヌツイスタは シフト演算を使って異なる位置のビットに影響を与えるようにした。 このアイディアにより、メルセンヌツイスタは、単一のビットのみならず、 32ビットで見てもk次均等分布するようになった。

サイコロを10回振って、全て1の目が現れた直後でも、 11回目に1の目が出る確率は6分の1であることは明らかだ。 メルセンヌツイスタは623次均等分布しているので 例えば0が連続して10回現れた後でも、 次にさらに0が現れる可能性がある。 メルセンヌツイスタは自然乱数の持つ性質を高レベルで実現している。

メルセンヌツイスタの周期は10進法で6千桁以上ある。 一般に「天文学的数字」といった場合でも、 せいぜい20桁くらいなので、この周期の長さは形容のしようがないほど長い。 また、GFSRでは、内部テーブルの初期化が非常に複雑であったのが、 メルセンヌツイスタでは全て0の場合を除けば、どのようにテーブルを初期化しても、 周期など理論的な性質に問題が起きない。また、生成速度もアルゴリズムが複雑なわりに、 掛け算や割り算など遅い命令を使わないので、意外と速い。

メルセンヌツイスタはTGFSRやTTGFSRと呼ばれていた時代から、 多くの研究者によって、テストされてきた。TGFSRの頃は問題が報告されていたが、 現在のメルセンヌツイスタに統計的な欠陥は見つかっていない。 筆者もいろいろテストしたが、問題は見つからなかった。

SFMT はメルセンヌツイスタの新しいものである。 SSE2 に対応した CPU で効率よく乱数を生成することができる。 また、旧メルセンヌツイスタでは、 ゼロを多く含む部分列が近くに集まりやすいという欠点があったが、SFMT では改善されている。 ただし、旧メルセンヌツイスタのこの欠点は、理論上の問題であり、 普通に種を与えて初期化していれば、統計上、全く問題は起きない。 なぜなら、非常に長い周期から見れば、極短い部分列での問題であり、 意図的に内部を操作しなければ、そのような部分列が現れることはまずない。

XorShift

高速な乱数に XorShift というのがある。

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)) ); 
} 

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

周期はメルセンヌツイスタほどはないが、実用上は十分といえる。 とりあえず、1億個の生成速度を XorShift を含めて Core 2 Duo E6600 + VS2005 で計ったところ、

 zsfmt(SSE2)   :  156ms
 SFMT(SSE2)    :  156ms
 zmtrand(SSE2) :  250ms
 zxor          :  265ms
 zmtrand       :  266ms
 xor128()      :  328ms
 zsfmt         :  375ms
 mt19937ar     :  625ms
 SFMT          : 1046ms
 rand()        : 1969ms
となった。SSE2 を使わない場合高速である。

論文の 2 ページ目に n が 32 と 64 の場合には T=(I+La)(I+Rb)では (a,b) が見つからなかったと書いてあるが、 n が 64 の場合 (7,9) は条件を満たしている。 また、同じく 2 ページ目の下にある 81 個の (a,b,c) のうち、|9,5,1| は |9,5,14| の間違いである。 また、3 ページ目の

(1) y^=y<<a; y^=y>>b; y^=y<<c;
(2) y^=y<<c; y^=y>>b; y^=y<<a;
(3) y^=y>>a; y^=y<<b; y^=y>>c;
(4) y^=y>>c; y^=y<<b; y^=y>>a;
(5) y^=y<<a; y^=y<<c; y^=y>>b;
(6) y^=y<<c; y^=y<<a; y^=y>>b;
(7) y^=y>>a; y^=y>>c; y^=y<<b;
(8) y^=y>>c; y^=y>>a; y^=y<<b;
の(5)と(6)は全く同じ結果となる。(7)と(8)も全く同じ結果となる。 さらに 4 ページ目の最初にあるプログラムの y=(y>>17); は y^=(y>>17); の間違いである。

http://www.iro.umontreal.ca/~lecuyer/myftp/papers/xorshift.pdf には、問題があると書かれている。

おまけ乱数 Xorshift のSSE2を使った高速化(整数版) ( http://www001.upp.so-net.ne.jp/isaku/xor.c.html )
乱数 Xorshift のSSE2を使った高速化(浮動小数点数版) ( http://www001.upp.so-net.ne.jp/isaku/dxor.c.html )
乱数 Xorshift の SSE2 を使わない高速化 ( http://www001.upp.so-net.ne.jp/isaku/fxor.c.html )