疑似乱数

疑似乱数を以前から使っているのですが、これはなかなか奥深い世界のようで、あまり考えないで作ると、どこかである種の規則性を持込んでしまいがちです。例えば、
　　　x := (x + a) mod n
といった数列は、a が素数である場合には、0 <= x < n の範囲を隈なく経めぐるものであることは確実ですが、隣接する、あるいは近接する数相互の関係には一定のくせがあるのが問題です。実際、こうした「乱数」によって点の位置を決め、プロットした結果は次のようになります。

より信頼できるとされるのは、「線形合同（linear congruential）法」と言われる次のものらしいです。
　　　x := (x * a + c) mod n

さらにシンプルな「乗算合同（multiplicative congruential）法」
　　　x := (x * a) mod n
というのもあり、「ニューメリカルレシピ（Numerical Recipes；Press,Teukolsky, Vettering, Flannery）」によると、「これ（乗算合同法）が「線形合同（linear congruential）法」と同程度に良いものでありうるという証拠は、理論的にも経験的にもある」とされているので、前回のページを通じて 1 <= x < n の範囲を全て埋める条件を確認できる乗算合同法だけについて検討しました。

n は素数（前々回の記事の手段で見つけられる）であると共に、a は最大周期のもの（前回の記事の手段で選別できる）であることが理想ですが、それは簡単に見つけられます。「ニューメリカルレシピ」では、n = 2147483647、a = 16807 を勧めていますが、ここでは、n = 161803398874991、a = 61803398874989 としてみました。上位桁と下位桁とを別けて利用する場合には、より信頼できるものになるのではないかなと思います。結果は、0 以上 1 未満、小数点以下 14 桁の実数として得られます。

C# と javascript で同じ結果になるルーチンを作りました。 (x * a) mod n という計算は一旦桁数が非常に大きくなるのが問題で、C# では decimal を使い、また、javascript では 10 進の桁ごとに、掛算と法を求めてシフトするやり方を採用しました。両方の結果を下図に示します。点の大きさが違って見えますが、詳細に見て頂くと、本質的に同じ結果と言えると思います（canvas による描画を追加しました --- 2025/1/1）。