Ｃ言語標準ライブラリの乱数rand( )は質に問題があり、禁止している学会もある。 他にも乱数には様々なアルゴリズムがあるが、多くのものが問題を持っている。 最も多くの人に使われている乱数であろう Visual Basic の Rnd の質は最低である。
乱数とは、本来サイコロを振って出る目から得られるような数を意味する。 このような乱数は予測不能なものである。 しかし、計算機を使って乱数を発生させた場合、 次に出る数は完全に決まっているので、予測不能とはいえない。 そこで、計算機で作り出される乱数を疑似乱数(PRNG)と呼び区別することがある。 ここでは、特にことわらない限り乱数とは疑似乱数のことを指すとする。
計算機でソフト的に乱数を発生させることの最大のメリットは、 再現性があることである。 初期状態が同じであれば、発生する乱数も全く同じものが得られる。 このことは、パラメータを変更して実験する場合、 乱数による揺らぎを完全に無くすのに便利である。 もちろん、サイコロを振って得られる乱数には再現性がない。
ほとんどの乱数生成プログラムは、種と呼ばれる数値を受け取って、 それに対応する乱数列を発生する。 実験のとき、種を別のものにすることによって、 簡単に異なる乱数列を得ることができる。
「乱数はもともとでたらめなので、いい加減な方法を使っても問題は起きない」 という考えがある。この考えは、危険である。 例えば、現在C言語の入門書などで広く紹介されている rand( ) という乱数は、 期待するほどでたらめではなく、ある意味で規則的な数列を返す。 従って、例えば rand( ) を使って実験を行ない、 統計処理した場合、結果の信頼性が低くなる。
このページの目的は、従来の乱数の問題点を指摘し、 すばらしい乱数、 「メルセンヌツイスタ」 の普及を促進することである。
また、乱数のことについて深く解説した杉田洋氏のページ モンテカルロ法，乱数，および疑似乱数も紹介しておく。
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の繰り返しになる。 すなわち、偶数と奇数が交互に生成される。 このことから、この乱数で下位ビットを乱数として使うのは危険であることがわかる。 また、ある乱数が得られたら、次に現れる乱数が１種類しかないという欠点も持つ。
多くの 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( )その１、その２に共通していることは、線形合同法を使っていることである。 すなわち、変数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( )その２ はX[n]の下位16ビットを捨て、 上位16ビットを右シフトしてそのうちの下位15ビットを返す。 従って最下位ビットでも周期は131072になるので、 「rand( ) その１」のような下位ビットを使った極端な問題は起きない。 しかし、乱数の返す数値の範囲が、0から32767までしかない。 言い換えれば、精度が15ビットしかない。 そのため、例えば、Box-Muller法で正規分布に従う乱数を発生させたり、 逆関数法で指数分布に従う乱数を発生させたりするとき、 端の方で理論的な分布と、生成された乱数の分布との違いが大きくなる。 下の表は全周期で最大値や最小値を調べた結果である。
正規分布 指数分布
最小値 最大値 最大値
rand()その１ -6.53σ 6.35σ 21.49λ
rand()その２ -4.56σ 4.56σ 10.40λ
rand( )の一番大きな欠点は、移植性がないことである。 rand( )その１の最大値は 2147483647 であるのに対し、 rand( ) その２は 32767 なので、異なるコンパイラでプログラムの移植するときに、 ソースコードを書き換えなくてはならない。 また、rand( )その２の間でもパラメータが違うものがあるので、 一方のコンパイラで起きたことを、他方のコンパイラで再現することが不可能である。
さらに、rand( )に共通する欠点として、違う種をsrand( )に与えても、 同じ乱数列が発生する場合があることも挙げられる。 例えば、srand(0)とsrand(0x80000000)で初期化された乱数列は、互いに同じものになる。
Excel などで使われる VBA の Rnd は簡略化して書けば以下のようになっている。 ※Excel2019 でも確認した。
static long x=327680;
float Rnd(void)
{
x=x*16598013+12820163&16777215;
return x*(1.0/16777216.0);
}
見てわかるように、24ビット線形合同法を使っている。 周期はたったの 1677万しかない。ちょっとした実験なら一周してしまううえ、精度も24ビットしかない。
もっと大きな問題は、乱数を初期化する Randomize にある。
これは与えられた浮動小数点数の種を16ビットの整数に変換して初期化する。
従って、最初に Randomize を行う場合、65536種類の系列しか得られない。
違う種を指定して実験しても、かなりの確率で、まったく同じ系列が返される危険がある。
また、下位 8 ビットは以前の下位 8 ビットがそのまま繰り越される。従って Randomize(1) として
計算を実行した後、やり直すために再び Randomize(1) として計算をしても、同じ結果を得られない。
参考までに、Rnd と Randomize と全く同じ動作をする MyRnd と MyRandomize を以下に示す。 (以前はＣ言語で作成した .dll と VBA で書いていたが、2021/1/30 に全て VBA に書き直した。)
自前のRnd ( http://isaku-wada.my.cooca.jp/rand/rnd.html )
Excel2010 以降のワークシート関数 RAND() はメルセンヌツイスタなので、これを使えば 品質の高い結果を得ることができる。
RAND 関数 - Office サポート
Range("A1").Value = Evaluate("RAND()")
のようにすると A1 セルに乱数を設定することができる。同じことを省略して、
[A1] = [RAND()]
のように書くこともできる。
ただし、種を指定する方法が無いので、得られる乱数列を再現することはできない。
Ｃ言語なら乱数のプログラムを書いて、対処できるが、
VBA の場合、整数のオーバーフローを認めていないので、種が指定できる性能のよい乱数を使いたいなら、
Ｃ言語などでダイナミック・リンク・ライブラリのプログラムを書いて、それを呼び出すしかない。
多くの乱数は 32 ビット整数演算をしているのでエクセルＶＢＡの Long 型では計算できない。
以下のプログラムは後述する Xorshift とメルセンヌツイスタを 31 ビットにしてＶＢＡで書いたものである。
ＶＢエディターを開いて、標準モジュールを追加し、そこに貼り付ければ使うことができる。 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で１億個発生させるのに Xorshift で約１分、メルセンヌツイスタで約５分かかった。
スピードを求めるなら乱数ライブラリ ( http://isaku-wada.my.cooca.jp/rand2.html ) のダイナミック・リンク・ライブラリ(libZD.iso)を使うとよいだろう。
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( )は２の４８乗の周期を持つ。 この周期の長さは rand( ) に比べて十分長い。 ただし、mrand48( )の最下位ビットの周期は 131072 しかないことに注意すること。 また、64ビットの掛け算を行っているため、 実行時間が rand( ) に比べかなり遅くなっている。
drand48( )の長所として、移植性の高さが挙げられる。 ほとんどの UNIX で drand48( ) が使用でき、全く同じ動作をする。 Visual C++で drand48( ) が使用できないのが残念に思われる。
rand( ) その1、rand( ) その２、drand48( ) は線形合同法で乱数を 発生させている。この方法に共通する欠点は、ある乱数が得られたら、 次に現れる乱数が限られてしまうことだ。とくに、rand( ) その１と、 drand48( ) は、次の乱数が完全に決まってしまう。 この欠点は乱数を組にして使うと問題になる。 例えば２次元ベクトル(x,y)を線形合同法で発生させると、 限られたベクトルしか現れない。
ちなみに、メルセンヌツイスタでは６２３次元のベクトルを発生させても、 全周期で、全ての組み合わせのベクトルが現れる。
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で、 全体の周期は約２の６２乗である。 random( )の優れた点は生成するときの実行時間が短いことである。 乱数の質が問題にならず、速度だけが問題になる時は、この乱数を使うとよい。
しかし、random( ) には統計的な欠陥が見つかっているので注意が必要である。 中心極限定理により、乱数を多数加えると正規分布に従うはずであるが、 random( ) を使った場合は正規分布から少し外れてしまう。
また、新しい srandom() の場合、種を０にした場合と１にした場合に、 全く同じ結果となる。この事は man 3 random、man 3 random_r と得られるマニュアルに 書かれておらず、ソースコードのコメントにも書かれていない。マルチスレッドで、 srandom_r( ) にスレッドの番号を渡すようなプログラムを組むような場合注意が必要である。
random( )のマニュアルには全てのビットがランダムで、 使用可能と書かれているが、以前の srandom( ) の初期化がまずいため、 下位ビットを使うのが好ましくないことがある。
例えば、以下のプログラムは、以前の random( ) を使って種を４つおきにして初期化し、 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( ) では４種類の全く同じ系列が順番に現れることがわかる。 この問題は現在の random( ) では解決されている。 しかし、この問題はlrand48( )で生じる。 従って、種を変えて100回実験をしても、得られる結果は4通りしかない。 次の表に、各乱数の最下位ビットの特性をまとめてみた。
周期 系列の種類の数
rand()その１ 2 2
rand()その２ 131072 131072
random() 4294967294 4
lrand48() 262144 4
mrand48() 131072 2
これを見てわかることは、最下位ビットを使って実験するとき、 以前の random( ) や lrand48( ) を使うよりも rand( )その２の方が、 ましな結果が得られるということである。
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); }
例えば、以下のＣ＃プログラムは、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 の乱数は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
となり、０から１０２３までの種で初期化したにもかかわらず、 最初の乱数は０.７近辺しか現われない。
また、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は各ビットが独立して生成されていた。メルセンヌツイスタは シフト演算を使って異なる位置のビットに影響を与えるようにした。 このアイディアにより、メルセンヌツイスタは、単一のビットのみならず、 ３２ビットで見てもk次均等分布するようになった。
サイコロを１０回振って、全て１の目が現れた直後でも、 １１回目に１の目が出る確率は６分の１であることは明らかだ。 メルセンヌツイスタは６２３次均等分布しているので 例えば０が連続して１０回現れた後でも、 次にさらに０が現れる可能性がある。 メルセンヌツイスタは自然乱数の持つ性質を高レベルで実現している。
メルセンヌツイスタの周期は10進法で６千桁以上ある。 一般に「天文学的数字」といった場合でも、 せいぜい20桁くらいなので、この周期の長さは形容のしようがないほど長い。 また、GFSRでは、内部テーブルの初期化が非常に複雑であったのが、 メルセンヌツイスタでは全て０の場合を除けば、どのようにテーブルを初期化しても、 周期など理論的な性質に問題が起きない。また、生成速度もアルゴリズムが複雑なわりに、 掛け算や割り算など遅い命令を使わないので、意外と速い。
メルセンヌツイスタはTGFSRやTTGFSRと呼ばれていた時代から、 多くの研究者によって、テストされてきた。TGFSRの頃は問題が報告されていたが、 現在のメルセンヌツイスタに統計的な欠陥は見つかっていない。 筆者もいろいろテストしたが、問題は見つからなかった。
SFMT はメルセンヌツイスタの新しいものである。 SSE2 に対応した CPU で効率よく乱数を生成することができる。 また、旧メルセンヌツイスタでは、 ゼロを多く含む部分列が近くに集まりやすいという欠点があったが、SFMT では改善されている。 ただし、旧メルセンヌツイスタのこの欠点は、理論上の問題であり、 普通に種を与えて初期化していれば、統計上、全く問題は起きない。 なぜなら、非常に長い周期から見れば、極短い部分列での問題であり、 意図的に内部を操作しなければ、そのような部分列が現れることはまずない。
高速な乱数に 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/
周期はメルセンヌツイスタほどはないが、実用上は十分といえる。
論文の 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 には、問題があると書かれている。