Fortranで作るXorshift(疑似乱数生成)

この記事は $2^{128}-1$ 周期のアルゴリズムonlyの解説のみになります。
# Xorshift
[メルセンヌ・ツイスタ(通称:MT)](https://ja.wikipedia.org/wiki/%E3%83%A1%E3%83%AB%E3%82%BB%E3%83%B3%E3%83%8C%E3%83%BB%E3%83%84%E3%82%A4%E3%82%B9%E3%82%BF)以外に高速で乱数を生成出来ないかなと考えて調べたところこんな記事を発見した。

https://ja.wikipedia.org/wiki/Xorshift

驚いたのは演算が排他的論理和とビットシフトのみであることだ。メルセンヌ・ツイスタより周期性は劣るが普通に使っていて問題ないだろうと考える。
ありがたいことに、WikipediaにC言語のソースコードが置いてあるのでそれをFortranに変えて実行してみよう思う。

以下ソースコード
```fortran:xorshift.f90
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値を設定できます。

以下ソースコード
```fortran:xorshift_.f90
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に毎回違う値を入れると出力される値も変わるので、まぁこんなもんでしょうか。

コンパイルして実行すると
```text:一回目
         251278488319
         648182000204
         475080158015
        3815321119129
      517440615441335
```
```text:二回目
         251278488287
         648182000236
         475080157983
        3815321184409
      517440615441303
```
どうでしょうか？値もちゃんと変わっていますね。

## 次は範囲指定をしたい
よくC言語で出回っているのは余りをとる方法ですね。
しかし、それだと余りを出しているだけなので偏りが出てしまいます。

じゃどうすんの？

はい、範囲以内になるまでループさせます。ですが範囲内になるまで結構時間がかかるのでここでは余りのみの解説とします。(実装の仕方がダメなのかそもそのこのプログラム自体がおかしいのか...)

0~nまで数値を10回表示させるソースコード
```fortran:xorshift_1.f90
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までの範囲まで出力させることができたので一応成功ですね。＾＾；
もっと範囲内でランダムに乱数を生成させたいときはこの記事を参考にするといいかもしれないです。
[なぜ乱数ジェネレーターを使用しているときにモジュロバイアスがあると人々が言うのですか？](https://www.web-dev-qa-db-ja.com/ja/c++/%e3%81%aa%e3%81%9c%e4%b9%b1%e6%95%b0%e3%82%b8%e3%82%a7%e3%83%8d%e3%83%ac%e3%83%bc%e3%82%bf%e3%83%bc%e3%82%92%e4%bd%bf%e7%94%a8%e3%81%97%e3%81%a6%e3%81%84%e3%82%8b%e3%81%a8%e3%81%8d%e3%81%ab%e3%83%a2%e3%82%b8%e3%83%a5%e3%83%ad%e3%83%90%e3%82%a4%e3%82%a2%e3%82%b9%e3%81%8c%e3%81%82%e3%82%8b%e3%81%a8%e4%ba%ba%e3%80%85%e3%81%8c%e8%a8%80%e3%81%86%e3%81%ae%e3%81%a7%e3%81%99%e3%81%8b%ef%bc%9f/1066683422/)

おわり