岡本安晴

Up

疑似乱数の生成

 

　周期が約である疑似乱数を生成するパッケージRngstreams (P. L’Ecuyer, et al. 2002)のアルゴリズムが、R. L. Eubank and A. Kupresanin (2012). Statistical Computing in C++ and R. CRC Pressのp. 129に説明されている。このアルゴリズムをVisual C++ 2008でプログラミングしてみた。計算効率よりアルゴリズムをわかりやすく忠実にという観点からプログラムコードを以下のように用意してみた。クラス型RNP2to191は、このウェブサイトの後半で説明しているサンプルプログラムの圧縮ファイルPlotRN.zipに含まれているヘッダーファイルRNP2to191.hにおいて宣言されているので、圧縮ファイル名PlotRN.zipをクリックしてダウンロード解凍すれば自由に利用できる。なお、ネイティブ版（Win32 アプリケーション）も用意しているが、これはこのウェブサイトの後半で扱っている。Linux版はウェブサイト＜http://y-okamoto-psy1949.la.coocan.jp/CppBldr/rnp2to191/＞に用意した。周期が約と長いので、分割して独立な複数の乱数系列として利用することができるが、これについてはこのウェブサイトの後の方で説明している。

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

        class RNP2to191 {

                public:

                                int ic, ib;

 

                                __int64 x1_0, x1_1, x1_2, x1_3;

                                __int64 x2_0, x2_1, x2_2, x2_3;

                                __int64 c2p32minus209, c2p32minus22853;

                                double inv2p32minus208, cforzeroZn;

 

                                inline __int64 DivModP( __int64 a, __int64 m ){

                                                __int64 v = a % m;

                                                if (v < 0) v += m;

                                                return v;

                                }

 

                                RNP2to191(){

                                        ic = ib = 0;

 

                                        c2p32minus209 = __int64(0x0000000100000000) - 209;

                                        c2p32minus22853 = __int64(0x0000000100000000) - 22853;

                                        inv2p32minus208 = 1.0 / (Math::Pow(2.0, 32.0) - 208.0);

                                        cforzeroZn = (Math::Pow(2.0, 32.0) - 209.0) / (Math::Pow(2.0, 32.0) - 208.0);

                                       

                                        x1_2 = 1;

                                        x1_1 = 2;

                                        x1_0 = 3;

                                        x2_2 = 4;

                                        x2_1 = 5;

                                        x2_0 = 6;

                                }

 

                                RNP2to191( int icycle, int ibias ){

                                        ic = icycle;

                                        ib = ibias;

 

                                        c2p32minus209 = __int64(0x0000000100000000) - 209;

                                        c2p32minus22853 = __int64(0x0000000100000000) - 22853;

                                        inv2p32minus208 = 1.0 / (Math::Pow(2.0, 32.0) - 208.0);

                                        cforzeroZn = (Math::Pow(2.0, 32.0) - 209.0) / (Math::Pow(2.0, 32.0) - 208.0);

                                       

                                        x1_2 = 1;

                                        x1_1 = 2;

                                        x1_0 = 3;

                                        x2_2 = 4;

                                        x2_1 = 5;

                                        x2_0 = 6;

 

                                        if (ib > 0)

                                                for (int i = 0; i < ib; i++) uni0();

                                }

・

・

・

                                inline void next(){

                                        x1_3 = x1_2;

                                        x1_2 = x1_1;

                                        x1_1 = x1_0;

                                        x2_3 = x2_2;

                                        x2_2 = x2_1;

                                        x2_1 = x2_0;

                                        x1_0 = DivModP( 1403580 * x1_2 - 810728 * x1_3, c2p32minus209 );

                                        x2_0 = DivModP( 527612 * x2_1 - 1370589 * x2_3, c2p32minus22853 );

               

                                        return;

                                }

・

・

・

                                double uni(){

                                        for (int i = 0; i < ic - 1; i++) next();

 

                                        next();

                                        __int64 zn = DivModP( x1_0 - x2_0, c2p32minus209 );

                                        if (zn > 0){

                                                        return zn * inv2p32minus208;

                                        }else{

                                                        return cforzeroZn;

                                        }

                                }

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

icに２以上の値が設定されると、一様乱数の関数uni( )はic個ずつ離れた次の値を返す。ibの値は、ic個ずつの値を返す時の最初の出発点となる乱数の系列の位置である。したがって、N個の独立した乱数の系列が必要なときは、icにNを設定して、ibにそれぞれ0からN−1までのN個の異なる値を与えてN個の乱数生成のオブジェクトを用意すればよい。例えば、3個の独立な乱数生成は

        RNP2to191 rn1(3, 0), rn2(3, 1), rn3(3, 2);

というように用意できる。

　上の乱数生成クラス型を用いたプログラム例を作成した。このプログラムPlotRN.exeをクリックしてダウンロード実行すると図１の画面になる。

図１

 

画面右上の「Uniform」ボタンをクリックすると一様乱数が

                                 for (int i = 0; i < 10000; i++){

                                         int xpos = int(x0 + (rn->uni() - 0.5) * 0.95 * ph);

                                         int ypos = int(y0 + (rn->uni() - 0.5) * 0.95 * ph);

                                         g->DrawEllipse( pn, xpos-2, ypos-2, 4, 4 );

                                 }

と呼び出されて図２のような点の分布が描かれる。

図２

 

　乱数生成のオブジェクトは、フォーム生成時に

                                 rn = new RNP2to191();

と生成されている。ポインタ型変数rnは、Form1のメンバとして

                                 RNP2to191 * rn;

と宣言されている。

「Normal」ボタンをクリックすると正規乱数が

                                 for (int i = 0; i < 5000; i++){

                                         double vx, vy;

                                         rn->normalPair( vx, vy );

                                         int xpos = int(x0 + 0.1 * vx * ph);

                                         int ypos = int(y0 + 0.1 * vy * ph);

                                         g->DrawEllipse( pn, xpos-2, ypos-2, 4, 4 );

                                 }

と呼び出されて図３のような点の布置が描かれる。

図３

 

プログラムのソースコードファイルなどは圧縮ファイルとして用意した。圧縮ファイル名PlotRN.zipをクリックしてダウンロードすることができる。ダウンロードした圧縮ファイルはマウスの右ボタンでクリックして表示されるメニュから「解凍」あるいは「展開」を選んで解凍することができる。解凍したプログラムファイルはVisual C++ 2008で開くことができる。ただし、このウェブサイトの後半で説明している乱数の系列において現在の位置から個後ろに飛ぶ関数Jump2toj(j)は、後半で説明しているプログラムのファイルCheckJump.zip（CLR版）、CheckJumpWin32.zip（Win32版）の方に含まれている。

 

　上の乱数の周期は約とされている。関数uni( )を回呼び出すのにどれくらいの時間が掛かるのかを調べるため以下のプログラムを作成した。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

#include "stdafx.h"

#include <Windows.h>

#include "RNP2to191.h"

 

using namespace System;

using namespace rnp2to191;

 

int main(array<System::String ^> ^args)

{

        RNP2to191 rn;

        int t0 = ::GetTickCount();

        for (int i = 0; i < 100000000; i++){

                rn.uni();

        }

        int t1 = ::GetTickCount();

        Console::WriteLine(((t1 - t0)/1000).ToString() + "s elapsed.");

        double vt = ((t1 - t0) * Math::Pow(2.0, 191) / 100000000.0) / double(1000.0 * 60.0 * 60.0 * 24.0);

        Console::WriteLine("It will take about " + vt.ToString("E1") + "days to call uni() 2**191 times.");

        Console::ReadLine();

    return 0;

}

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

Dell XPS8300; Core i7 CPUでのデバッグモードでの実行結果を図４に示す。

図４

 

約日掛かるという予測である。系列内の乱数を全て呼び出し尽くすことは現実問題として不可能である。後ろの方の系列の乱数は使われることがないということになる。後ろの方の乱数を呼び出すために途中の乱数の呼び出しを飛び越えて後ろの方に飛ぶ関数を用意した。このウェブサイトの後の方で説明している。

「デバッグなし」で実行すると図５の結果になる。

図５

 

上の実行時間計時プログラムTimeRNP2to191.exeはクリックしてダウンロード実行することができる。プログラム・ソースコード・ファイルなどを圧縮したものも圧縮ファイル名TimeRNP2to191.zipをクリックしてダウンロードできる。ダウンロードした圧縮ファイルはファイル名TimeRNP2to191.zipをマウスの右ボタンでクリックして表示されるメニュから「解凍」あるいは「展開」を選んで解凍することができる。解凍したプログラムファイルはVisual C++ 2008で開くことができる。

 

　ネイティブ版（win32 コンソールアプリケーション）も用意した。上のサンプルコード（CLR版）では、ドットネットのコンポーネントを利用しているが、これをVisual C++ 2008のwin32対応のものに書き改めた。乱数生成のヘッダーファイルおよびソースコードファイルはMathクラスのメンバー関数をC++の標準のものに変更した。メインのソースコードファイルも以下のように変更した。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

#include "stdafx.h"

#include <iostream>

#include <iomanip>

#include "RNP2to191.h"                                                  //      乱数用ヘッダーファイル

#include "finetimer.h"                                                  //      計時用ヘッダーファイル

 

using namespace std;

using namespace rnp2to191;                                              //      乱数用名前空間

using namespace mytimer;                                                //      計時用名前空間

 

 

int _tmain(int argc, _TCHAR* argv[])

{

        RNP2to191 rn;                                                   //      乱数生成オブジェクト

        MyTimer mytimer;                                                //      計時オブジェクト

        int t0 = mytimer.timems();                                      //      ミリ秒単位の計時

        for (int i = 0; i < 100000000i64; i++){

                rn.uni();

        }

        int t1 = mytimer.timems();

        cout << (t1 - t0)/1000 << "s elapsed." << endl;

        double vt = ((t1 - t0) * pow(2.0, 191) / 100000000.0) / double(1000.0 * 60.0 * 60.0 * 24.0);

        cout << "It will take about " << setprecision(5) << vt << "days to call uni() 2**191 times." << endl;

        cout << "Press Enter-key.";

        cin.get();

       

        return 0;

}

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

ネイティブ版をデバッグモードで実行すると図６のようになった。

図６

 

CLR版より時間が掛かっている。Release版を作成して「デバッグなし」で実行すると図７の結果となった。

図７

 

ドットネット版の「デバッグなし」で実行より少し速くなっている。

ネイティブ版（win32アプリケーション版）のプログラムTimeRNP2to191Win32.exeも用意したので、プログラム名TimeRNP2to191Win32.exeをクリックしてダウンロード実行することができる。プログラム・ソースコード・ファイルなどをまとめたものも圧縮ファイルTimeRNP2to191Win32.zipとして用意した。圧縮ファイル名TimeRNP2to191Win32.zipをクリックしてダウンロード解凍すれば、Visual C++ 2008で開くことができる。この圧縮ファイルには、win32用の乱数生成ヘッダーファイルおよびソースコードファイルも含まれているが、上のCLR版に対応する部分を以下に示しておく。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

        class RNP2to191 {

                public:

                                int ic, ib;

 

                                __int64 x1_0, x1_1, x1_2, x1_3;

                                __int64 x2_0, x2_1, x2_2, x2_3;

                                __int64 c2p32minus209, c2p32minus22853;

                                double inv2p32minus208, cforzeroZn;

 

                                inline __int64 DivModP( __int64 a, __int64 m ){

                                                __int64 v = a % m;

                                                if (v < 0) v += m;

                                                return v;

                                }

 

                                RNP2to191(){

                                        ic = ib = 0;

 

                                        c2p32minus209 = __int64(0x0000000100000000) - 209;

                                        c2p32minus22853 = __int64(0x0000000100000000) - 22853;

                                        inv2p32minus208 = 1.0 / (pow(2.0, 32.0) - 208.0);

                                        cforzeroZn = (pow(2.0, 32.0) - 209.0) / (pow(2.0, 32.0) - 208.0);

                                       

                                        x1_2 = 1;

                                        x1_1 = 2;

                                        x1_0 = 3;

                                        x2_2 = 4;

                                        x2_1 = 5;

                                        x2_0 = 6;

                                }

 

                                RNP2to191( int icycle, int ibias ){

                                        ic = icycle;

                                        ib = ibias;

 

                                        c2p32minus209 = __int64(0x0000000100000000) - 209;

                                        c2p32minus22853 = __int64(0x0000000100000000) - 22853;

                                        inv2p32minus208 = 1.0 / (pow(2.0, 32.0) - 208.0);

                                        cforzeroZn = (pow(2.0, 32.0) - 209.0) / (pow(2.0, 32.0) - 208.0);

                                       

                                        x1_2 = 1;

                                        x1_1 = 2;

                                        x1_0 = 3;

                                        x2_2 = 4;

                                        x2_1 = 5;

                                        x2_0 = 6;

 

                                        if (ib > 0)

                                                for (int i = 0; i < ib; i++) uni0();

                                }

 

                                inline void next(){

                                        x1_3 = x1_2;

                                        x1_2 = x1_1;

                                        x1_1 = x1_0;

                                        x2_3 = x2_2;

                                        x2_2 = x2_1;

                                        x2_1 = x2_0;

                                        x1_0 = DivModP( 1403580 * x1_2 - 810728 * x1_3, c2p32minus209 );

                                        x2_0 = DivModP( 527612 * x2_1 - 1370589 * x2_3, c2p32minus22853 );

               

                                        return;

                                }

 

                                double uni(){

                                        for (int i = 0; i < ic - 1; i++) next();

 

                                        next();

                                        __int64 zn = DivModP( x1_0 - x2_0, c2p32minus209 );

                                        if (zn > 0){

                                                        return zn * inv2p32minus208;

                                        }else{

                                                        return cforzeroZn;

                                        }

                                }

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

乱数呼び出しの跳び越し

 

　乱数生成の基になっている数列は次式

によって関係づけられている。ここで、の二乗を回繰り返すとの乗が得られるので、個後の乱数が得られることになる。このAの二乗を繰り返して得られるによって個後に跳ぶ関数として

                void Jump2toj( int j );

を用意した。この関数Jump2tojの使用例を以下に示す。　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

        RNP2to191 rn, rn1, rn2, rn3;

        rn1.Jump2toj(2);

        rn2.Jump2toj(2); rn2.Jump2toj(2);

        rn3.Jump2toj(4);

        for (int i = 0; i < 20; i++){

                Console::WriteLine(i.ToString() + "   " + rn.uni().ToString() +

                                                        "   " + rn1.uni().ToString() + "   " + rn2.uni().ToString() +

                                                        "   " + rn3.uni().ToString() );

        }

        Console::WriteLine();

        Console::WriteLine("Press Enter-key.");

        Console::ReadLine();

    return 0;

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

実行結果は以下のようになる。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

0   0.00100949784041744   0.466827597259577   0.759799247383662   0.344027651603741

1   0.595003783879985   0.378907797116046   0.811096096110527   0.217561267840849

2   0.357834537613574   0.00693863431067111   0.0774230598714195   0.768123863211312

3   0.222340826701115   0.994035933343571   0.941320921945095   0.970742524115938

4   0.466827597259577   0.759799247383662   0.590508772252552   0.0503500142304234

5   0.378907797116046   0.811096096110527   0.665767405759455   0.719490179012985

6   0.00693863431067111   0.0774230598714195   0.367942868622978   0.739075256448158

7   0.994035933343571   0.941320921945095   0.232095737307312   0.347864491249391

8   0.759799247383662   0.590508772252552   0.344027651603741   0.768540035434144

9   0.811096096110527   0.665767405759455   0.217561267840849   0.304796551213992

10   0.0774230598714195   0.367942868622978   0.768123863211312   0.0142935689010337

11   0.941320921945095   0.232095737307312   0.970742524115938   0.167940764439218

12   0.590508772252552   0.344027651603741   0.0503500142304234   0.242746087836844

13   0.665767405759455   0.217561267840849   0.719490179012985   0.509156073188517

14   0.367942868622978   0.768123863211312   0.739075256448158   0.906678396879953

15   0.232095737307312   0.970742524115938   0.347864491249391   0.609260818857292

16   0.344027651603741   0.0503500142304234   0.768540035434144   0.726640496203029

17   0.217561267840849   0.719490179012985   0.304796551213992   0.173988858282027

18   0.768123863211312   0.739075256448158   0.0142935689010337   0.512383043434395

19   0.970742524115938   0.347864491249391   0.167940764439218   0.0384586013386466

 

Press Enter-key.

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

関数Jump2tojによる乱数の系列の飛越が、色を付けた文字を対応させることによって確認することができる。なお、Jump2toj(j)をk回呼び出すと、乱数の現在の位置から個後ろに飛ぶ。

上のプログラム（CLR版）をヘッダーファイルなどと一緒に圧縮ファイルとしてまとめたものCheckJump.zipも用意した。圧縮ファイル名CheckJump.zipをクリックしてダウンロードしたものはVisual C++ 2008で開くことができる。

Win32版は圧縮ファイルCheckJumpWin32.zipにまとめた。圧縮ファイル名CheckJumpWin32.zipをクリックしてダウンロードすれば解凍してVisual C++ 2008で開くことができる。

十分に離れたところに飛ぶことにより、実質的に独立な乱数系列を複数得ることができる。

関数uni()を回呼ぶのに必要な時間と個後の位置に飛ぶJump2toj(100)の実行時間を計るプログラムを以下のように作成した。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

#include "stdafx.h"

#include <Windows.h>

#include "RNP2to191.h"

 

using namespace System;

using namespace rnp2to191;

 

int main(array<System::String ^> ^args)

{

        RNP2to191 rn;

        int t0 = ::GetTickCount();

        for (int i = 0; i < 100000000; i++){

                rn.uni();

        }

        int t1 = ::GetTickCount();

        Console::WriteLine(((t1 - t0)/1000).ToString() + "s elapsed.");

        double vt = ((t1 - t0) * Math::Pow(2.0, 100) / 100000000.0) / double(1000.0 * 60.0 * 60.0 * 24.0);

        Console::WriteLine("It will take about " + vt.ToString("E1") + "days to call uni() 2**100 times.");

        t0 = ::GetTickCount();

        rn.Jump2toj(100);

        t1 = ::GetTickCount();

        Console::WriteLine("It took " + (t1 - t0).ToString() +"msec to execute Jump2toj(100).");

        Console::WriteLine();

        Console::WriteLine("Press Enter-key.");

        Console::ReadLine();

    return 0;

}

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

上のプログラムの実行結果例を図８に示す。

図８

 

関数uni()の回の呼び出しにはオーダーの日数が掛かるとの予想であるが、Jump2toj(100)の実行はGetTickCountの精度以下の時間である。上のプログラムは圧縮ファイルTimeJump.zipにまとめた。

また、Jump2tojのLinux版は別のウェブサイト＜http://y-okamoto-psy1949.la.coocan.jp/CppBldr/rnp2to191/＞で扱っている。

 

 

なお、Visual C++プログラミングの入門書として＜岡本安晴「大学生のための心理学VC++プログラミング入門」勁草書房＞を上梓した。

 

参考文献

L’Ecuyer, P., Simard, R., Chen, E. J. and Kelton, W. D. (2002). An object-oriented random-number package with many long streams and substreams. Operations Research, 50, 1073-1075.

 

Up