BASIC で乱数

BASIC で乱数

〜十進 BASIC で疑似乱数を生成〜

2022-10-10 作成福島

乱数とは

乱数とは、でたらめに現れる数のことです。
何かを決めるとき、自分では無い他の誰かに決めてほしいことが多々あります。

トランプをかき混ぜる時は誰かの意思が働かないほうが望ましいし、じゃんけんでは相手の手が分かっているとつまらなくなります。
ギャンブルにおいては、予め答えが分かっていたら勝負にならない上に、場合によっては罪に問われることもあります。

また、煙や水の動きなど、多すぎる不確定な数値をすべて人間の手によって決定するのは現実的ではありません。

そんな時に便利なのが「乱数」です。
自分の意志とは関係なく、でたらめに現れる数が答えを決めてくれます。

コンピュータにおける乱数

コンピュータにおける乱数は、実は乱数ではありません。
正しくは「疑似乱数」と呼び、一見して乱数に見える数字の羅列を計算で生成して、本物の乱数の代用にしています。
このため、アルゴリズムと初期値さえ同じであれば、同じ結果を得ることができます。

本物の乱数は物理法則に従うので、同じ結果を得ることはできません。
よく、「空中に放り投げたコインが着地した面の向き」や「サイコロを振った時の出目(でめ)」として使われます。
コインやサイコロを投げた時に全く同じ環境を再現することが物理的にできないことは理解できるでしょう。

環境が変わっても同じ乱数列を

最近のプログラミング言語には、ほぼすべてに高機能な乱数生成器が内蔵されており、内蔵でなくても付録として付いていたりします。
また、Linux の様に高機能な OS には乱数生成器が内蔵されており、気軽に利用できるようになっています。
アプリケーションを作るときはプログラムからこれらを使用すれば良いので、通常は乱数生成器そのものを作成する必要はありません。

十進 BASIC では「N = INT(RND * 10)」の記述で N に 0～9 の乱数を求めることができます。

便利な機能ですが、欠点もあります。
プログラミング言語や OS に内蔵されている乱数発生器は、それぞれ独自のアルゴリズムによって実装されているため、環境が異なると同じ乱数列を発生させることができません。
また、アルゴリズムが同一でも実数を返す形式の乱数発生器だと、小数の丸め誤差により値が同一になりません。

具体的には、乱数をキーにした暗号器で暗号化した文書を他の言語で復号化できなかったり、
シミュレーションの結果が、異なる OS で同じ再現にならないことになります。

ここでは乱数発生器を作成し、異なる環境でも同じ乱数列を生成させることにより、同じ実行結果を得られるようにします。

疑似乱数のアルゴリズム

ここでは、Xorshift という手法の中から 128bit 版の xorshift128() を実装します。
これは周期が 2⁹⁶-1 になる (らしい^*1) のです。
^*1すみません、ちゃんと証明していません。

疑似乱数の周期とは、生成した数字と同じパターンが次に現れるまでの回数のことを指します。
例えば、
1, 3, 2 , 1, 3, 2 , 1, 3, 2 , 1, 3, 2 , ···
であれば 3、
2, 13, 3, 11, 7 , 2, 13, 3, 11, 7 , 2, 13, 3, 11, 7 , ···
なら 5 です。

計算で疑似的に乱数を生成するので、どうしても周期が現れてしまいます。
2⁹⁶-1 は想像しにくいので 10 進数に直すと、

95 ≦ log₂X ＜ 96
95 ≦ log₁₀X / log₁₀2 ＜ 96
95 ≦ log₁₀X / 0.3010 ＜ 96
28.595 ≦ log₁₀X ＜ 28.896
2⁹⁶-1 = 10^{28.595～28.895}-1

となり、10 進数に直しても想像しにくいほど長い周期をもつということが分かりました。(10²⁸ は 1 穣(じょう))
(1 秒間に 100 億回の乱数を生成すると仮定して約 2,500 億年必要。えーと、ビッグバン説が正しいとして、宇宙の誕生って 138 億年前だっけ)

xorshift128 には制約があって、内部的に符号なし 32bit の整数演算が出来なければなりません。
十進 BASIC は符号なし 32bit の整数演算を直接することはできませんが 15 桁までの実数が扱えるので、
これを使って符号なし 32bit の整数演算を実現します。(2³²⁺¹ は 10¹⁰ に収まる)

余談: VBA の場合、整数は符号あり 32bit の演算しかできませんが、複数の Byte 型変数を連結することにより、符号なし 32bit の実装が可能です。(四則演算は面倒だけど、論理演算とシフトだけなら結構カンタン)

乱数発生器の実装

アルゴリズムは xorshift128() の、そのまんまパクリです。
(符号なし 32bit 演算が前提なので、直接に負数や実数を扱うことは出来ません)

DIM xorshift_state(4)
DATA 123456789, 0, 0, 0     ! すべてが 0 にならないように初期化する。
MAT READ xorshift_state
DECLARE EXTERNAL SUB Xorshift128

FOR i = 1 TO 100            ! 10²⁹ 個も確認してられないので、とりあえず 100 個。
   CALL Xorshift128(xorshift_state)
   LET n = xorshift_state(1)
   PRINT n  ! 乱数を表示する。(実際に使うときは MOD(n, x) であまりを求める^*2)
NEXT i
END


! 乱数発生器
EXTERNAL SUB Xorshift128(state())
LET t = MOD(state(4), 2^32)     ! fool proof^*3
LET s = MOD(state(1), 2^32)     ! fool proof^*3
LET state(4) = state(3)
LET state(3) = state(2)
LET state(2) = s
LET t = MOD(bitXOR(t, t * 2^11), 2^32)              ! t ^= t << 11
LET t = bitXOR(t, INT(t / 2^8))                     ! t ^= t >> 8
LET state(1) = bitXOR(bitXOR(t, s), INT(s / 2^19))  ! t ^ s ^ (s >> 19)
END SUB

^*2例えば 0～9 の値を求めるには「MOD(n, 10)」と記述する。
^*3この 2 箇所の MOD(n, 2^32) は、人間が初期値として n < 2³² を暗算するのが苦手のため記述している。本来は不要。

実行結果
このプログラムと Python 版の実行結果が同一になることを確認してください。

おまけ

乱数を使いやすくするために、公差が 1 の等差数列をかき混ぜる副プログラム DustShuffle を掲載します。(Dust: Richard Durstenfeld)
ダステンフェルドのシャッフルは、本来は要素を限定しないのですが、この副プログラムでは等差数列としています。
実行結果の数列を添え字として使用することで、本来の目的を果たせるように作ってあります。
DustShuffle の中から呼び出している副プログラム Xorshift128 は上記と同じものです。
DIM xorshift_state(4)
DATA 123456789, 0, 0, 0     ! すべてが 0 にならないように初期化する。
MAT READ xorshift_state
DECLARE EXTERNAL SUB Xorshift128
DECLARE EXTERNAL SUB DustShuffle

DIM idx(8)
CALL DustShuffle(idx, SIZE(idx), xorshift_state)

FOR i = 1 TO 8
   PRINT idx(i);
NEXT i
PRINT

END


! 乱数発生器
EXTERNAL SUB Xorshift128(state())
LET t = MOD(state(4), 2^32)     ! fool proof
LET s = MOD(state(1), 2^32)     ! fool proof
LET state(4) = state(3)
LET state(3) = state(2)
LET state(2) = s
LET t = MOD(bitXOR(t, t * 2^11), 2^32)              ! t ^= t << 11
LET t = bitXOR(t, INT(t / 2^8))                     ! t ^= t >> 8
LET state(1) = bitXOR(bitXOR(t, s), INT(s / 2^19))  ! t ^ s ^ (s >> 19)
END SUB


! ダステンフェルドのシャッフル
EXTERNAL SUB DustShuffle(arr(), siz, state())
FOR i = 1 TO siz
   LET arr(i) = i
NEXT i
FOR i = siz TO 1 STEP -1
   CALL Xorshift128(state)
   LET p = MOD(state(1), i) + 1
   LET n = arr(p)
   LET arr(p) = arr(i)
   LET arr(i) = n
NEXT i
END SUB
本当は引数 siz を省略したいところだけど、十進 BASIC は配列の長さを可変にできないためこうしています。

実行結果
 6  4  8  5  1  3  2  7 