疑似乱数２

**デフォルトの名無しさん** · 2007/10/17(水) 22:34:59

擬似乱数発生器について語ろうか。その２

前スレ
擬似乱数
http://pc11.2ch.net/test/read.cgi/tech/1146071975/

関連スレ
【危険】とんでもプログラム告発スレッド【悪質】
http://pc11.2ch.net/test/read.cgi/tech/1191860116/

SIMD-oriented Fast Mersenne Twister (SFMT):
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index.html

**デフォルトの名無しさん** · 2007/11/01(木) 01:20:36

>>29
ttp://www.jstatsoft.org/v08/i14/
ttp://www.iro.umontreal.ca/~lecuyer/myftp/papers/xorshift.pdf

**デフォルトの名無しさん** · 2007/11/02(金) 15:50:53

ド素人なので突っ込みどころがあったら御教授くださいorz
rand()で生成した乱数をバイナリ形式で出力して
NIST SP800-22やdiehardテストに突っ込みたいんですけど・・・

↓のままだとdiehardどころかSP800-22にも引っかかってしまいます
#include <stdio.h>
#include <stdlib.h> // rand, srand使用
#include <time.h> // time使用
#define size 12000000 //bit列の長さ定義

↓

**デフォルトの名無しさん** · 2007/11/02(金) 15:52:35

main()
{
FILE *outputfile; // 出力ストリーム
unsigned long int m,i=0;
char bit[size];
srand((unsigned) time(NULL)); // time関数からシードをセット
outputfile = fopen("bit.dat", "w"); // ファイルを書き込み用にオープン
if (outputfile == NULL) { // オープンに失敗した場合
printf("cannot open\n"); // エラーメッセージを出して
exit(1); // 異常終了
}
for(i=0; i<size-1; i++){ // 乱数を発生させ剰余を計算
m=rand()%0xffff;
bit[i]=m;
}
fwrite(&bit,size,1,outputfile); //バイナリの書き込み
fclose(outputfile); // ファイルをクローズ
return 0;
}

↓公式ページ
NIST ttp://csrc.nist.gov/groups/ST/toolkit/rng/documentation_software.html
diehard ttp://stat.fsu.edu/pub/diehard/

**デフォルトの名無しさん** · 2007/11/03(土) 03:37:01

>>31
rand()は乱数としての性質があまりよろしくないという
当たり前のことを文句言われてもどうかと。

**デフォルトの名無しさん** · 2007/11/03(土) 12:45:15

>>32

> rand()%0xffff

一般論として、
剰余を取っても良質な疑似乱数列となっていることが保証されている
乱数生成系でない場合は、剰余ではなく商を取って一定の範囲に
切り詰めるべき。

歴史的に、システムのデフォルトの乱数生成系は品質の悪いものが
多いので、良い疑似乱数が必要なら、マニュアル等で良質な
乱数生成系を使っていることが確認できなければ、デフォルトの
生成系は使うべきでない。

**デフォルトの名無しさん** · 2007/11/03(土) 17:26:32

そんなのどうでもいいじゃん
だからまあ普通はラッパーかけて開発中はrand使って
あとで精度ほしくなったら差し替えるように作るわけだが

31 · 2007/11/03(土) 23:00:03

>>33-35
早速の返答ありがとうございます。
色々とマヌケな勘違いをしていましたorz

rand()は「暗号用乱数に相応しくない乱数」の一例として示すために用いました。
剰余をとったのは、良質なRNGならば下位bitの乱数性もまた良質だろう
ということを逆説的に示すためです。

NISTに同梱されている線形合同法のRNGはテストをパスし、
逆に上記のプログラムで生成した乱数列は全く話にならなかったので、
バイナリへの変換の仕方がマズかったのかと思っていました。

NISTは線形合同法くらいなら通ってしまうと聞いたので
てっきりrand()でもイケるものだと…orz

**デフォルトの名無しさん** · 2007/11/04(日) 05:32:17

>>36
実行はできてるんだね。
オレの環境だと>>31-32のコードは動かないんだけど。
# 原因はスタックに置くには配列サイズがでかすぎることみたいだけど。

乱数の質に絡みそうな部分というと型の選択がまずいんじゃね？
charじゃなくてunsigned charの配列にしたらどう？
それと剰余が違うんじゃね？
rand()%0xffff → rand()&0xffff でしょ。
あとループの回数も一回少ないと思う。

**デフォルトの名無しさん** · 2007/11/05(月) 02:14:18

RAND__MAXと0xffffってどっちが大きいか知ってる？

**デフォルトの名無しさん** · 2007/11/10(土) 21:14:12

RAND_MAXの値は実装にもよるかと。

そんなことより、rand()%0xffff → rand()&0xffff　　←こいつに誰か突っ込まんのか。

**デフォルトの名無しさん** · 2007/11/10(土) 21:27:10

うん

**デフォルトの名無しさん** · 2007/11/10(土) 21:39:02

そんなことよりxorshiftの話しよーぜ

**デフォルトの名無しさん** · 2007/11/11(日) 19:11:45

unsigned long xor32(){
static unsigned long y=2463534242;
yˆ=(y<<13);y=(y>>17);return (yˆ=(y<<5));}

unsigned long long xor64(){
static unsigned long long x=88172645463325252LL;
xˆ=(x<<13);xˆ=(x>>7);return(xˆ=(x<<17));}

unsigned long xor96(){
static unsigned long x=123456789,y=362436069,z=521288629;unsigned long t;
t=(xˆ(x<<20))ˆ(yˆ(y>>11))ˆ(zˆ(z<<27))ˆ(wˆ(w>>6));x=y;y=z;z=w;return(w=t);}

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

unsigned long xor160(){
static unsigned long x=123456789,y=362436069,z=521288629,w=88675123,v=5783321;unsigned long t;
t=(xˆ(x>>7));x=y;y=z;z=w;w=v;return v=(vˆ(v>>6))ˆ(tˆ(t>>13));}

unsigned long xor192(){
static unsigned long x=123456789,y=362436069,z=521288629,w=88675123,v=5783321,d=6615241;unsigned long t;
t=(xˆ(x>>2));x=y;y=z;z=w;w=v;v=(vˆ(v<<4))ˆ(tˆ(t<<1));return(d+=362437)+v;}
こんな具合か

42 · 2007/11/11(日) 19:16:33

コピペ丸張りで色々とおかしくなってることに気付く。
後悔はしていない。

**デフォルトの名無しさん** · 2007/11/11(日) 20:03:25

>>40
意味とランダム度が全然違うだろ。

**デフォルトの名無しさん** · 2007/11/11(日) 20:05:36

>>42
これってsrand()に相当する関数は自作しなきゃいかんよね？

**デフォルトの名無しさん** · 2007/11/13(火) 04:38:48

>>45
srand_xor(int hoge)
{
extern static int x;
x=hoge;
}
これでいいんとちゃう？ん？staticはexternできなかったっけ？

**デフォルトの名無しさん** · 2007/11/13(火) 08:50:49

>>46
アホかい。

**デフォルトの名無しさん** · 2007/11/13(火) 17:51:10

分からなかったらとりあえずクラス化→private属性のクラス内変数で（ｒｙ

**デフォルトの名無しさん** · 2007/11/16(金) 23:27:42

あるシードを使うと、得られる乱数を繋げたバイナリが
偶然たまたま迷走Mind.mp3 として読めるような
乱数ジェネレータを作ったらどうだろうか。

**デフォルトの名無しさん** · 2007/11/16(金) 23:44:09

確率的にありえない

**デフォルトの名無しさん** · 2007/11/17(土) 00:09:54

意図的に作るんなら偶然じゃないな

**デフォルトの名無しさん** · 2007/11/17(土) 00:27:25

PSG

**デフォルトの名無しさん** · 2007/11/17(土) 05:19:51

>>42 には間違いがあるので注意

**デフォルトの名無しさん** · 2007/11/17(土) 05:24:13

訂正したコードを示した方が有用なレスになると思うよ。

**デフォルトの名無しさん** · 2008/01/02(水) 19:12:39

>>42 の

unsigned long xor32(){
static unsigned long y=2463534242;
y?=(y<<13);y=(y>>17);return (y?=(y<<5));}

の y=(y>>17) は y^=(y>>17) とすべき

**デフォルトの名無しさん** · 2008/01/06(日) 13:04:45

XORってライブラリ化してないの？

**ヽ・´∀｀・,,）っ━━━━━━┓** · 2008/01/07(月) 00:17:55

めぼしいのはBoost.Randomにある

**デフォルトの名無しさん** · 2008/01/07(月) 12:16:11

団子あげ

**デフォルトの名無しさん** · 2008/01/13(日) 16:23:59

MTの別実装を作ってみました。
ttp://mt-lite.sourceforge.net/index.html.ja

基本のアルゴリズムはいじらず実装だけを作ったもので、
メモリ節約・リアルタイム性・スレッドセーフなどに一通り配慮したつもりです。

mt19937ar・SFMT・WELL・Xorshiftなどとの比較評価はこのあたりに:
ttp://mt-lite.sourceforge.net/doc/ja/evaluation.summary.html

**デフォルトの名無しさん** · 2008/01/14(月) 01:05:03

乙

**ヽ・´∀｀・,,）っ━━━━━━┓** · 2008/01/14(月) 01:28:22

SFMTのスレッドセーフ実装なら京大の人のページにあったぞ。SSE2だけだけど。

**デフォルトの名無しさん** · 2008/01/14(月) 01:48:05

XorShiftをseed使えるようにする場合、どうすればいいかな？＜何度もブン回す以外で

**デフォルトの名無しさん** · 2008/01/14(月) 04:04:29

http://web.comlab.ox.ac.uk/oucl/work/richard.brent/random.html

**デフォルトの名無しさん** · 2008/01/14(月) 17:24:18

>>63
thx

**デフォルトの名無しさん** · 2008/01/14(月) 21:52:52

/* xor128()をSSE2で実装してみました。XorShift は SIMD に向かない事がわかりました。*/
#include <stdio.h>
#include <time.h>
#include <emmintrin.h>
unsigned long xor128(void){
static unsigned long x=123456789UL,y=362436069UL,z=521288629UL,w=88675123UL;
unsigned long t;t=(x^(x<<11));x=y;y=z;z=w;return(w=(w^(w>>19))^(t^(t>>8)));}

unsigned long xor128sse2(void){
static union{unsigned long v[4];__m128i m;}u={123456789UL,362436069UL,521288629UL,88675123UL};
__m128i x=u.m, w, t, r; r=_mm_slli_epi32(x,11);t=_mm_xor_si128(x,r); w=_mm_srli_si128(x,12);
x=_mm_srli_si128(x,4); r=_mm_srli_epi32(w,19);w=_mm_xor_si128(w,r); r=_mm_srli_epi32(t,8);
t=_mm_xor_si128(t,r); w=_mm_xor_si128(w,t); r=_mm_slli_si128(w,12);x=_mm_or_si128(x,r);
_mm_store_si128(&u.m,x);return(u.v[3]);}

void main(void){long i,c;
for(i=0;i<16;i++)printf("%8lx %8lx\n",xor128(),xor128sse2());
c=clock();for(i=0;i<50000000L;i++)xor128(); printf("xor128(): %4lu msec\n",clock()-c);
c=clock();for(i=0;i<50000000L;i++)xor128sse2();printf("xor128sse2():%4lu msec\n",clock()-c);}

**デフォルトの名無しさん** · 2008/01/14(月) 22:04:20

32 行書けるんだからもうちょっと綺麗におながいします

**デフォルトの名無しさん** · 2008/01/14(月) 23:52:54

上のほうでrand()をダイハードに食わせるとかいってた人がいたのでオチを教えておきます。
p=1.0頻発します。テストの種類によっては部分がすべてp=1.0　あわせてp=1.0みたいな。

xor128（）って、例の論文に載ってるコードを組み込む場合、ライセンスはPD扱いですか？

**デフォルトの名無しさん** · 2008/01/15(火) 00:53:58

>>65
こりゃ無茶だろ．こういう使い方は．．．
別の種の乱数を４個同時に生成するのが正解

**デフォルトの名無しさん** · 2008/01/15(火) 16:03:59

>>62
xor128内部で使われてる変数x,y,z,wの中のxに、デフォルトの数値(123456789)の変わりに
与えられたseed値を代入するようにしてみた。
一応ちゃんと動作してるっぽいが統計的に良いのか悪いのかはよく分からん

**デフォルトの名無しさん** · 2008/01/15(火) 18:39:47

>>69
初期値ってあの値じゃなくても良かったのね。
乱数の質がどの程度変わるのかが気になるが

**デフォルトの名無しさん** · 2008/01/15(火) 19:28:29

やっつけでtimeの値入れたり

**デフォルトの名無しさん** · 2008/01/15(火) 19:41:19

SEED値を入れるならxよりwの方がいい気がする

**デフォルトの名無しさん** · 2008/01/15(火) 20:24:10

/* >>68 さん。的確なアドバイスありがとうございます。>>68 さんの方法で書いてみました。*/
/* 種で初期化する関数付き。XorShift の SSE2 による高速化はあきらめるしかないのか！ */
#include <stdio.h>
#include <time.h>
#include <emmintrin.h>
unsigned long xor128(void) {
static unsigned long x=123456789UL,y=362436069UL,z=521288629UL,w=88675123UL;
unsigned long t;t=(x^(x<<11));x=y;y=z;z=w;return(w=(w^(w>>19))^(t^(t>>8)));}

int idx=4; union { unsigned long v[4]; __m128i m; }
x={123456789,123456789,123456789,123456789},y={362436069,362436069,362436069,362436069},
z={521288629,521288629,521288629,521288629},w={ 88675123, 88675123, 88675123, 88675123};

void sxor128x4(unsigned long s)
{int i; for (idx=4,i=0;i<16;i++) x.v[i]=s=1812433253*(s^(s>>30))+i; }

unsigned long xor128x4(void) { __m128i t,r,s;
if (idx==4) { r=_mm_slli_epi32(x.m,11); t=_mm_xor_si128(x.m,r); x.m=y.m; y.m=z.m;
s=z.m=w.m; r=_mm_srli_epi32(s,19); s=_mm_xor_si128(s,r); r=_mm_srli_epi32(t,8);
t=_mm_xor_si128(t,r); s=_mm_xor_si128(s,t); w.m=s; idx=0; }
return(w.v[idx++]); }

void main(void){ long i,c;
for (i=0;i<9;i++) printf("%8lx %8lx %8lx %8lx %8lx\n",
xor128(),xor128x4(),xor128x4(),xor128x4(),xor128x4() );
for (sxor128x4(1),i=0;i<9;i++)
for (printf("\n"),c=0;c<4;c++) printf("%8lx ",xor128x4() );
for (c=clock(),i=0;i<100000000L;i++) xor128();
printf("\n(xor128():%lu msec), ",clock()-c );
for (c=clock(),i=0;i<100000000L;i++) xor128x4();
printf("(xor128x4():%lu msec)\n",clock()-c ); }

**デフォルトの名無しさん** · 2008/01/15(火) 21:20:55

>>73
その種の使い方、MTと全く同じだな

**デフォルトの名無しさん** · 2008/01/17(木) 22:51:40

XorShiftの正しいシードの設定方法を教えてくれ！

**デフォルトの名無しさん** · 2008/01/26(土) 18:55:32

論文の６ページ目に xor128 の seed set は全て０の場合を除いて、どの
ように設定してもよいと書いてある。しかし、極端な設定をすると不自然
な部分列が出力される。そこで初期化の方法はＭＴを参考にした。ｉを０
からではなく１から始めた理由は０を種にしたときに最初の２つの出力が
似た値にならないようにするためである。配列を使っているが、スピード
はオリジナルとほとんど変わらない。コンパイラによっては、こちらの方
が速くなる場合もある。デフォルトの初期ベクトルは rand,srand になら
って種を１としたときと同じである。

unsigned long MyVec128[4]=
{ 1812433254UL, 3713160357UL, 3109174145UL, 64984499UL };

void MySeed128(unsigned long s)
{
int i;
for (i=1;i<=4;i++) MyVec128[i-1]=s=1812433253UL*(s^(s>>30))+i;
}

unsigned long MyXor128(void)
{
unsigned long *a=MyVec128, t=a[0], w=a[3];
a[0]=a[1]; a[1]=a[2]; a[2]=w;
t^=t<<11; t^=t>>8; w^=w>>19; w^=t; a[3]=w; return w;
}

**デフォルトの名無しさん** · 2008/01/29(火) 19:35:39

なるほど、どんな初期値であっても、2^128-1のどこかの部分列になってるのか

**デフォルトの名無しさん** · 2008/01/29(火) 23:17:42

>>61について誰か詳しく頼む

**デフォルトの名無しさん** · 2008/01/31(木) 21:25:58

ん？　何が分からなかった？

**デフォルトの名無しさん** · 2008/01/31(木) 21:44:19

遅いんだよ
もういい

**デフォルトの名無しさん** · 2008/01/31(木) 23:24:35

はやっ

**デフォルトの名無しさん** · 2008/01/31(木) 23:53:55

こんな過疎スレで遅い言われても

**デフォルトの名無しさん** · 2008/01/31(木) 23:59:56

自分がいくら張り付いてるからって、ね。

78 · 2008/02/02(土) 06:55:22

>>61のソースの所在教えてくれ
パクる

**ヽ・´∀｀・,,）っ━━━━━━┓** · 2008/02/02(土) 07:20:15

dSFMTだったごめん。
http://nct.brain.riken.jp/~takekawa/softwares.xhtml

あとこんなのもあった。
http://charles.karney.info/random/

**デフォルトの名無しさん** · 2008/02/03(日) 13:34:07

thx

**デフォルトの名無しさん** · 2008/02/03(日) 13:46:53

ダンゴさんはSSEネタには造詣が深いな。

**デフォルトの名無しさん** · 2008/02/03(日) 16:46:36

じゃあ上げとこ

**デフォルトの名無しさん** · 2008/02/22(金) 20:01:38

シミュ板、数学板に乱数スレ発見。
乱数
http://science6.2ch.net/test/read.cgi/sim/1100375806/
乱数の生成方法を必死になって考えるスレ　
http://science6.2ch.net/test/read.cgi/math/1203153071/

**デフォルトの名無しさん** · 2008/03/23(日) 01:46:24

あげ

**デフォルトの名無しさん** · 2008/03/27(木) 23:56:06

XORSHIFTって、その出力の剰余を取った場合も
良質なランダム性があるって言われているのでしょうか？

**デフォルトの名無しさん** · 2008/03/28(金) 00:05:34

XORSHIFTは良質ではないと思います

**デフォルトの名無しさん** · 2008/03/28(金) 00:06:25

工工エエエｴｴ（´Д｀）ｴｴエエエ工工

**デフォルトの名無しさん** · 2008/03/28(金) 00:11:48

MTが一番良質と思います　XORは速度は速いです　用途によっては十分です

**デフォルトの名無しさん** · 2008/03/28(金) 00:20:58

>>91です。
XORSHIFTで剰余を取った時の結果を載せてる
WEBページを見た気がするのですが、結果が
今一だったような記憶が・・・

自分でも試してみようと思っているのですが、
理論的には説明出来そうにないし。

**デフォルトの名無しさん** · 2008/03/28(金) 03:02:22

>>95
そのWEBページのURL教えてくれ

95 · 2008/03/28(金) 23:26:21

>>96
個人のブログですが、確かここだったと思います。
ttp://d.hatena.ne.jp/Isoparametric/20070830/1188462129

今、改めて読んでみると、よく意味分からん・・・

**デフォルトの名無しさん** · 2008/03/28(金) 23:55:04

たった100回回しただけなら偏りも出るだろう

**デフォルトの名無しさん** · 2008/04/15(火) 04:27:31

はぐれめたるは捕まえられますか？

**デフォルトの名無しさん** · 2008/04/17(木) 03:21:30

>>99
むり

**デフォルトの名無しさん** · 2008/04/20(日) 18:00:50

単純に何の補正も無い相手に対して「捕まえられる」確率はどのくらいか。
そしてはぐれメタルが逃げない確率はどのくらいか。
その両方のデータをくれ。

**デフォルトの名無しさん** · 2008/04/20(日) 19:05:05

乱数とどういう関係があるんだ

**デフォルトの名無しさん** · 2008/04/20(日) 23:27:07

え？はぐれメタルの行動を決めている乱数テーブルの解析の話じゃないの！？

**デフォルトの名無しさん** · 2008/04/21(月) 00:54:50

ここは「擬似乱数生成アルゴリズム」に関する議論のスレッドだ

**ヽ・´∀｀・,,）っ━━━━━━┓** · 2008/04/21(月) 13:20:43

PS2のドラクエ５では壁の衝突回数をカウントして乱数列の初期化してるようだ。
オラクルベリーの教会から一度も壁に衝突せずにカジノのスロットに直行すると
全く同じパターンが出てくるとか。

おそらくはぐれメタル捕獲にも何らかのパターンがあるかと

**デフォルトの名無しさん** · 2008/04/23(水) 08:14:14

全然ちげーよ、ばかが

**デフォルトの名無しさん** · 2008/04/23(水) 12:16:04

「ヽ・´∀｀・,,）っ━━━━━━┓」は放置推奨クソコテ

**デフォルトの名無しさん** · 2008/05/01(木) 23:36:04

英語版MTを読んでいたらこんなのを見つけた。
ttp://en.wikipedia.org/wiki/Multiply-with-carry
で、cで実装してみた。
あんまり使えないと思うけど。

/* Complementary Multiply With Carry */
void initialize( unsigned long );
unsigned long cmwc( void );
enum { A = 3636507990, B = 0xffffffff, R = 1024 };
unsigned long seed[ R ], index, c;

void initialize( unsigned long n )
{
unsigned long i;
for (i=0; i<R; i++) {
seed[ i ] = n;
n = n * 1103515245 + 12345;
}
index = c = 0;
}
unsigned long cmwc( void )
{
unsigned long x;
unsigned long long t;
t = (unsigned long long)seed[ index ] * A + c;
x = B - (unsigned long)( t & B );
c = (unsigned long)( t >> 32 );
seed[ index ] = x;
if ( ++index >= R ) index = 0;
return x;
}

**デフォルトの名無しさん** · 2008/05/03(土) 21:01:47

>>105
http://www.astr.tohoku.ac.jp/~ajiki/cheap_restaurant/GAME/SAGA2/saga2rand.html
サガ2の乱数発生器

DQ5もソース見れば多分判ると思うよ

**デフォルトの名無しさん** · 2008/05/18(日) 22:07:45

DebianのOpenSSLライブラリに予測可能な乱数の生成を行う脆弱性 - ITmedia
ttp://www.itmedia.co.jp/enterprise/articles/0805/15/news023.html
↓
OpenSSLの脆弱性突くブルートフォース攻撃発生、簡単に暗号解読の恐れ - ITmedia
ttp://www.itmedia.co.jp/news/articles/0805/16/news019.html

Debian JPプロジェクトのアナウンス
ttp://www.debian.or.jp/blog/openssl_package_and_its_vulnerability.html

**デフォルトの名無しさん** · 2008/05/19(月) 09:17:00

まとめ乙

**デフォルトの名無しさん** · 2008/08/06(水) 01:14:02

保守

**デフォルトの名無しさん** · 2008/08/09(土) 20:00:17

M系列という単語を見かけたので、理解してないまま保守リンク
http://pc11.2ch.net/test/read.cgi/tech/1158367586/751-

**デフォルトの名無しさん** · 2008/08/11(月) 14:05:33

VBAで、dll使わずに高機能な乱数を実装する方法はないですか？
XorShiftを使うには、長桁計算しかないですか？

**デフォルトの名無しさん** · 2008/08/13(水) 17:47:46

' Xorshift を VBA で書く場合、Double を使うと比較的シンプルになる
' VBA の演算子 Xor と関数 Fix は 31 bit までしか扱えないので
' 53 bit まで扱える関数 uFix、uXor を用意した。
Option Explicit
Public x, y, z, w As Double

Function uFix(ByVal x As Double) As Double
Const Base = 2# ^ 22
Dim y As Long
y = Int(x / Base): uFix = y * Base + Int(x - y * Base)
End Function

Function uXor(ByVal i As Double, ByVal j As Double) As Double
Const Base = 2# ^ 22
Dim u, v As Long
u = Int(i / Base): v = Int(j / Base): i = i - u * Base: j = j - v * Base
uXor = (u Xor v) * Base + (Int(i) Xor Int(j))
End Function

Function xor128() As Double
Dim t As Double
t = x * 2# ^ 11: t = t - Int(t / 2# ^ 32) * 2# ^ 32
t = uXor(x, t): x = y: y = z: z = w: t = uXor(t, uFix(t / 2# ^ 8))
w = uXor(uXor(w, uFix(w / 2# ^ 19)), t): xor128 = w
End Function

Sub Test()
Dim i As Long
x = 123456789#: y = 362436069#: z = 521288629#: w = 88675123#
For i = 1 To 2000: Cells(i, 1) = xor128: Next
End Sub

**115** · 2008/08/13(水) 18:47:55

' 書き込んだ後、uFix は必要がないことに気づきました。
' こっちの方が高速です。
Option Explicit
Public x, y, z, w As Double

Function uXor(ByVal i As Double, ByVal j As Double) As Double
Const Base = 2# ^ 30
Dim u, v As Long
u = Int(i / Base): v = Int(j / Base): i = i - u * Base: j = j - v * Base
uXor = (u Xor v) * Base + (Int(i) Xor Int(j))
End Function

Function xor128() As Double
Dim t As Double
t = x * 2# ^ 11: t = t - Int(t / 2# ^ 32) * 2# ^ 32
t = uXor(x, t): x = y: y = z: z = w: t = uXor(t, Int(t / 2# ^ 8))
w = uXor(uXor(w, Int(w / 2# ^ 19)), t): xor128 = w
End Function

Sub Test()
Dim i As Long
x = 123456789#: y = 362436069#: z = 521288629#: w = 88675123#
For i = 1 To 2000: Cells(i, 1) = xor128: Next
End Sub

**デフォルトの名無しさん** · 2008/08/18(月) 18:18:26

>>116
１１４です、ありがとうございます。32bitのXorShiftですね。

で、周期の途中から乱数を取り出して使いたいのですが、
x,y,z,wの初期値を乱数で設定するのに、デフォルトのrnd()を使ったら無意味ですよね。

**デフォルトの名無しさん** · 2008/08/19(火) 00:23:17

>>75あたりか

**115** · 2008/08/22(金) 21:40:15

'種を使って初期化できます。116よりさらに２倍程度高速になりました
Private Const p32 = 2# ^ 32, p31 = 2# ^ 31, p21 = 2# ^ 21, p11 = 2# ^ 11
Private Const m53 = 2# ^ -53, m32 = 2# ^ -32, m30 = 2# ^ -30, m19 = 2# ^ -19, m11 = 2# ^ -11, m8 = 2# ^ -8
Private x, y, z, w As Double
Private f As Boolean
Private Function uXor(ByVal x As Double, ByVal y As Double) As Double
Dim u, v As Long
If x >= p31 Then u = x - p32 Else u = x
If y >= p31 Then v = y - p32 Else v = y
u = u Xor v
If u < 0 Then uXor = u + p32 Else uXor = u
End Function
Private Function XSub(ByVal x As Double, ByVal i As Long) As Double
Dim s As Variant
s = CDec(1812433253) * CDec(uXor(x, Int(x * m30))) + CDec(i): s = s - CDec(Int(s * m32)) * CDec(p32): XSub = s
End Function
Public Sub InitXor(ByVal s As Long) ' s を種にして乱数を初期化する
f = True: x = XSub(s, 1): y = XSub(x, 2): z = XSub(y, 3): w = XSub(z, 4)
End Sub
Private Function NextXor() As Double
Dim t As Double
If f = 0 Then InitXor (1)
t = x * p11: t = t - Int(t * m32) * p32: t = uXor(x, t): x = y: y = z: z = w: t = uXor(t, Int(t * m8)): w = uXor(uXor(w, Int(w * m19)), t): NextXor = w
End Function
Public Function NextUnif() As Double ' 0 以上 1 未満の乱数を返す
Dim x, y As Double
x = NextXor * m11: y = NextXor: NextUnif = (y * p21 + Int(x)) * m53
End Function

**デフォルトの名無しさん** · 2008/08/22(金) 23:22:43

xorshiftって周期とか分布の保証ってあるの？
もしないならM系列の方がましだと思うんだけど。
ググれば長周期の多項式も見つかるし。

**デフォルトの名無しさん** · 2008/08/23(土) 01:50:28

劣化Ｍ系列と理解してるが

**デフォルトの名無しさん** · 2008/08/23(土) 02:50:13

早くて軽くてrandよりマシ

**デフォルトの名無しさん** · 2008/08/25(月) 16:31:11

>>119
excelVBAじゃ初期化するための種が結局16bitしかないんだから
Xorshiftの種も16bit種類しかなくならね？

**119** · 2008/08/27(水) 03:32:41

種 s は Long なので32bit種類あるが、負の場合は試してないので実質31bit種類。

**デフォルトの名無しさん** · 2008/08/27(水) 10:43:09

VBAのRandomizeが16bit種類
それを使って初期化後rnd()で取得できる乱数も16bit種類
XorShiftの種にそれらを使うのなら結局系列は16bit種類

本当の最初の最初に32bitの乱数をセットするときにどうすればいいのか。
32bitのランダムな位置からXorShiftを開始するのに、そのうちの16bit分しか設定できない。

**デフォルトの名無しさん** · 2008/08/27(水) 10:56:10

よく分からんが服を買いにいく服がない状態か

**119** · 2008/08/27(水) 14:14:51

>>125

確かにそうだ。困ったもんだ。

**デフォルトの名無しさん** · 2008/08/27(水) 14:18:21

>>125
GetTickCountとかCryptGenRandomとかWin32APIに頼ればいいじゃない。

**デフォルトの名無しさん** · 2008/08/29(金) 12:44:39

>>128
とりあえず
Private Declare Function GetTickCount Lib "kernel32" () As Long
後に
　　GetTickCountを使って取得するようにしたよ。