この記事は $2^{128}-1$ 周期のアルゴリズムonlyの解説のみになります。
Xorshift
メルセンヌ・ツイスタ(通称:MT)以外に高速で乱数を生成出来ないかなと考えて調べたところこんな記事を発見した。
驚いたのは演算が排他的論理和とビットシフトのみであることだ。メルセンヌ・ツイスタより周期性は劣るが普通に使っていて問題ないだろうと考える。
ありがたいことに、WikipediaにC言語のソースコードが置いてあるのでそれをFortranに変えて実行してみよう思う。
以下ソースコード
program xorshift
implicit none
integer(8) i
do i = 1, 5 ! 5回ループ
print *, xor128()
end do
contains
integer(8) function xor128() ! Xorshiftのアルゴリズム
implicit none
integer(8), save :: x=123456789, y=362436069, z=521288629, w=88675123
integer(8) t
t = xor(x, lshift(x, 11))
x = y; y = z; z = w;
w = xor(xor(w, rshift(w, 19)), xor(t, rshift(t, 8)))
xor128 = w
end function xor128
end program xorshift
このソースコードは、Xorshiftアルゴリズムで生成された乱数を5回表示するというものです。
試しにコンパイルして実行すると下のようになりました。
252977563114
646616338854
476657867818
294684809458
517442357777055
おお乱数だ！じゃ、2回目やってみるか！
やってみますと同じ値が出力されました。(;´д｀)トホホ
毎回違う乱数を生成するには？
結論はseed値を毎回変えることです。どうすりゃええんや？
考えに考えた結果、ここではtime関数を使って毎回違う疑似乱数を生成していこうと思います。まずtime関数は1970年1月1日からの経過秒数を返す関数で、これを利用すると簡単に毎回違うseed値を設定できます。
以下ソースコード
program xorshift
implicit none
integer(8) i, seed
seed = time() ! time()がtime関数で、経過時間をseedに代入
do i = 1, 5
print *, xor128(seed) ! xor128にseedを渡す、そして出力された値を表示
end do
contains
integer(8) function xor128(w)
implicit none
integer(8), intent(inout) :: w ! ここが変わると値も変わる
integer(8), save :: x=123456789, y=362436069, z=521288629
integer(8) t
t = xor(x, lshift(x, 11))
x = y; y = z; z = w;
w = xor(xor(w, rshift(w, 19)), xor(t, rshift(t, 8)))
xor128 = w
end function xor128
end program xorshift
コードがざっと2行増えたところでしょう。
変数wに毎回違う値を入れると出力される値も変わるので、まぁこんなもんでしょうか。
コンパイルして実行すると
251278488319
648182000204
475080158015
3815321119129
517440615441335
251278488287
648182000236
475080157983
3815321184409
517440615441303
どうでしょうか？値もちゃんと変わっていますね。
次は範囲指定をしたい
よくC言語で出回っているのは余りをとる方法ですね。
しかし、それだと余りを出しているだけなので偏りが出てしまいます。
じゃどうすんの？
はい、範囲以内になるまでループさせます。ですが範囲内になるまで結構時間がかかるのでここでは余りのみの解説とします。(実装の仕方がダメなのかそもそのこのプログラム自体がおかしいのか...)
0~nまで数値を10回表示させるソースコード
program xorshift
implicit none
integer(8) i, seed, n
seed = time()
read *, n ! nまでの数値を入力
do i = 1, 10
print *, mod(xor128(seed), n + 1) ! ここが肝
end do
以下同じなので省略
実行結果(10を入力したとき)
10
0
9
3
8
8
9
5
7
10
4
0~nまでの範囲まで出力させることができたので一応成功ですね。＾＾；
もっと範囲内でランダムに乱数を生成させたいときはこの記事を参考にするといいかもしれないです。
なぜ乱数ジェネレーターを使用しているときにモジュロバイアスがあると人々が言うのですか？
おわり