counter

Visual Basicサンプル

合同法による一様乱数の生成


乱数の利用方法には大きく分けて2種類あるものと思われます。ひとつは暗号に応用するために、生成される数列を予想することができない、という意味の乱数です。もうひとつは、数列は予想できてよい、数列が予め明示的に与えられていてさえよい、しかし統計的観点から十分にランダムな数列であること、少なくとも、統計的な制限範囲がよく研究され明らかな数列であることを要求する乱数です。厳密な定量的な統計処理を必要とする分野では後者のタイプの乱数が望まれます。


ここでは、特に後者のタイプの応用に適する乱数生成法を、Visual Basicで記述したソースコードで紹介します。また、ランダムさには問題がある(上述のとおりどれぐらい問題があるかは良く知られている)ものの、すばやく乱数を生成できる方法も説明します。


乗算合同法による一様乱数

VB.NET版 download randmulcon.NET.zip


Ran0、Ran1、及びRan2と3種類のクラスがありますが、いずれも


NUMERICAL RECIPES in C (日本語版)
William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery
技術評論社
ISBN4-87408-560-1


に掲載されている関数を移植したもので、31bit整数の乱数を内部的に生成し、これを基に0〜1のDouble型の一様乱数を生成するメソッドを備えています。(Double型で乱数が取り出されますが、内部的には高々31bitの精度しかありません。)



応用にあたっての留意・注意事項

大方、出典のNumerical Recipesに書かれている内容の繰り返しになりますが、Visual Basicの特徴も配慮して、応用の際に留意あるいは注意して頂きたいことを掲げておきます。擬似乱数については、様々な誤解の流布や、個々の開発者の思い込みによって、誤用が多い分野ですので、参考にして頂ければ幸いです。



インラインで実現できる32bit線形合同法(.Net 2.0のみ可能)

32bit符号無しの整数型(UInteger型あるいはUint32型)を用いて、かつ、整数のオーバーフローチェックを無効にした場合、以下のコードで、統計的な振る舞いは良くないですが非常に高速に乱数列を得られます。

Dim idum As UInteger ' = 何か初期値を代入する。0でもよい。 '何らかの処理のループ idum = 1664525UI * idum + 1013904223UI 'idumの値を乱数として利用する。 'idumの値はループ内の他のところで変更しないこと。 'ループの終わり


上述のidumの値を基にSingle型の浮動小数点数の乱数を得たい場合には、可読性には欠けますが、以下のような構造体のトリックを使って、やはりインラインで変換することができます。まずは、構造体SingleUnionの宣言を示します。

'構造体のメンバー変数としてSingle型とUInteger型(共に同じ '4Byte長)をメモリ上の同じ位置を共有させることができる。 'CやC++で言うところの共用体(union)と同等である。 'この機能は.Net Frameworkの以下の名前空間で提供されている。 Imports System.Runtime.InteropServices '属性を用いてメンバー変数のメモリ上の位置を細かく指定できる。 '詳しくは、.Net Frameworkのドキュメントを参照のこと。 <StructLayout(LayoutKind.Explicit)> _ Public Structure SingleUnion <FieldOffset(0)> Public S As Single <FieldOffset(0)> Public UI As UInteger End Structure


この構造体を用いて、UInteger型の乱数値をSingle型に変換します。UInteger型でビット操作を行いSingle型の内部表現に変換しますが、この操作がそのままSingle型のメンバー変数の値に反映されます。なお、IEEEに準拠するSingle型にのみ適用可能ですのでご注意ください。

'変換操作を行うSingleUnionの変数を宣言する。 Dim r As SingleUnion '上述のidumを算出する線形合同法の式 idum = 1664525UI * idum + 1013904223UI 'UInteger型のメンバー変数を9bit右シフトすることにより、idumの '上位23bitがSingle型の仮数部に置かれる。さらに、指数部に1以上、 '2未満の数値を示す値をセットする。 r.UI = &H3F800000UI Or (idum >> 9) 'この時点で、r.Sは[1,2)の範囲の乱数になるので、1だけ減じて '[0,1)に収まるようにする。 r.S -= 1.0! 'あとは、r.Sを利用する。


最初に掲げたサンプルのソースコードでコンパイルオプションを選択することにより、ここに示したものと同じ(ただし、31bitの整数を8bit右シフトする)方法を、最初に掲げたRan0、Ran1、及びRan2についても適用することができます。


改定日 内 容
2007/06/07 乗算合同法のサンプルを初掲載
2007/06/08 コンパイルオプションによって、内部で保持されている31bit整数の乱数の上位23bitを、直接Single型変数の仮数部にセットできる処理も選べるように、サンプルを変更。ファイル名も変更。
2007/06/19 乱数をDouble型で返すように変更。ただし、精度は高々31bitのままである。初期値を省略して初期化を行った場合、コンパイルオプションによっては、内臓時計の値に基づいて初期値を算出するように変更。繰り返し同じ系列の乱数を生成できるように、初期値を返すプロパティーを追加。