岡本安晴

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に説明されている。このアルゴリズムをLinux (fedora 17/64ビット)のC++でプログラミングしてみた。Fedora 17/64ビットでは、long int型は６４ビットの整数型である（long int型の有効ビット数を調べるプログラムはこのウェブサイトの後半で説明している）。計算効率よりアルゴリズムをわかりやすく忠実にという観点からプログラムコードを以下のように用意してみた。なお、乱数の周期が長いので分割して独立な複数の乱数系列として利用することができるが、これについてはこのウェブサイトの後の方で説明している。

なお、３２ビットOSのLinux用のプログラムは、このウェブサイトの最後に用意した。

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

#include <cmath>

 

namespace rnp2to191 {

 

              class RNP2to191 {

                            public:

                                                        int ic, ib;

 

                                                        long int x1_0, x1_1, x1_2, x1_3;

                                                        long int x2_0, x2_1, x2_2, x2_3;

                                                        long int c2p32minus209, c2p32minus22853;

                                                        double inv2p32minus208, cforzeroZn;

 

                                                        inline long int DivModP( long int a, long int m ){

                                                                                    long int v = a % m;

                                                                                    if (v < 0) v += m;

                                                                                    return v;

                                                        }

 

                                                        RNP2to191(){

                                                                      ic = ib = 0;

                                                                      c2p32minus209 = 0x0000000100000000L - 209;

                                                                      c2p32minus22853 = 0x0000000100000000L - 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 = 0x0000000100000000L - 209;

                                                                      c2p32minus22853 = 0x0000000100000000L - 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();

                                                                      long int zn = DivModP( x1_0 - x2_0, c2p32minus209 );

                                                                      if (zn > 0){

                                                                                                  return zn * inv2p32minus208;

                                                                      }else{

                                                                                                  return cforzeroZn;

                                                                      }

                                                        }

 

                                                       

             //    Standard Normal Distribution

              //    Rejection polar method for normal variate

                            double normal();

 

 

           //    Normal Distribution of mean m  and standard deviation s

                            double normalMS(double m, double s);

 

 

           //    A Pair of Independent Random Numbers of

              //      the Standard Normal Distribution

                            void  normalPair( double & n1, double & n2 );

 

 

              };

 

}

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

 

乱数のオブジェクト生成時に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);

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

上の乱数の利用例を以下に示す。

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

#include <iostream>

#include "RNP2to191.h"                                                                                       //            random number generator

 

using namespace std;

using namespace rnp2to191;                                                                               //            random number generator

 

 

int main()

{

              cout << "The program started." << endl;

              RNP2to191 rn, rn1(2, 0), rn2(2, 1), rnn;                                               //            random number generators

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

                                          cout << rn.uni() << "  " << rn1.uni() << "  " << rn2.uni() <<

                                                         "  " << rnn.normal() << endl;

              }

             

              return 0;

}

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

 

上のプログラムの実行結果は図１のようになる。

図１

 

ソースコードファイル、ヘッダーファイルなどを圧縮ファイルCheckRN.zipにまとめて用意した。圧縮ファイル名CheckRN.zipをクリックしてダウンロードすれば、自由に利用できる。ただし、このウェブサイトの後半で説明している乱数の系列において現在の位置から個後ろに飛ぶ関数Jump2toj(j)は、後半で説明しているプログラムのファイルCheckJumpLinux.zip（64ビットOS用）、CheckJumpLinux32.zip（32ビットOS用）の方に含まれている。

 

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

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

#include <iostream>

#include <iomanip>

#include <ctime>

#include "RNP2to191.h"                                                                                       //            random number generator

 

using namespace std;

using namespace rnp2to191;                                                                               //            random number generator

 

 

int main()

{

              cout << "The program started." << endl;

              RNP2to191 rn;                                                                                                      //            random number generator

              clock_t vt0 = clock();

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

                                          rn.uni();

              }

              clock_t vt1 = clock();

              cout << double(vt1 - vt0) / double(CLOCKS_PER_SEC) << "s elapsed." << endl;

              double vt = ((vt1 - vt0) * pow(2.0, 191) / 100000000.0) / (double(CLOCKS_PER_SEC) * 60.0 * 60.0 * 24.0);

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

             

              return 0;

}

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

実行結果を図２に示す。

図２

 

約日掛かるという予測である。この結果は、VM VirtualBox上にインストールされたfedora 17/64ビットでのものである。約日掛かるということは、乱数の系列の前の方以外は、実際上は呼び出されないということになる。したがって、系列の後ろの方に跳ぶことにより、１つの系列を分割して、それぞれの分割された部分系列を独立な乱数生成系列として用いることができる。系列の後ろの方に跳ぶ方法については、このウェブサイトの後半で説明している。

上の実行時間計時プログラムは、ソースコード・ファイル、ヘッダーファイルをまとめて圧縮したものCheckTime.zipを用意した。圧縮ファイル名CheckTime.zipをクリックしてダウンロード解凍して自由に利用することができる。

 

 

long int型の有効ビット数

 

　long int型の有効ビット数は、例えば、次のプログラムで調べることができる。

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

#include <iostream>

 

using namespace std;

 

int main(){

                            long int a = 1;

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

                                          a *= 2;

                                          cout << (i+2) << ":  " << a << endl;

                            }

                            return 0;

}

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

上のプログラムをfedora17/64ビットで実行すると図３のような結果を得る。

図３

 

６５ビット目でオーバーフローが起こり、０になっていることがわかる。すなわち、有効桁数は６４ビットである。

同じプログラムをfedora16/32ビットで実行すると図４の結果となる。

図４

 

３３ビット目でオーバーフローが起こり、０になっているので、有効桁数が３２ビットであることがわかる。

 

乱数生成の跳び越し

 

　乱数の生成において個後ろに飛ぶ関数Jump2toj(j)を用意した。使用プログラム例と実行結果を以下に示す。実行結果では、跳び越し先と対応する数値を同じ色で示した。なお、Jump2toj(j)をk回呼び出すと、乱数の系列の現在の位置から個後ろに飛ぶ。

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

[yasuharu@localhost CheckJumpLinux]$ cat CheckJumpLinux.cpp

 

#include <iostream>

#include "RNP2to191.h"

 

using namespace std;

using namespace rnp2to191;

 

int main()

{

              RNP2to191 rn, rn1, rn2, rn3;

              rn1.Jump2toj(2);

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

              rn3.Jump2toj(4);

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

                            cout << i << "   " << rn.uni() << "   " << rn1.uni()

                                          << "   " << rn2.uni() << "   " << rn3.uni() << endl;

              }

              cout << endl;

              return 0;

}

 

[yasuharu@localhost CheckJumpLinux]$ g++ -o check CheckJumpLinux.cpp RNP2to191.cpp

[yasuharu@localhost CheckJumpLinux]$ ./check

0   0.0010095   0.466828   0.759799   0.344028

1   0.595004   0.378908   0.811096   0.217561

2   0.357835   0.00693863   0.0774231   0.768124

3   0.222341   0.994036   0.941321   0.970743

4   0.466828   0.759799   0.590509   0.05035

5   0.378908   0.811096   0.665767   0.71949

6   0.00693863   0.0774231   0.367943   0.739075

7   0.994036   0.941321   0.232096   0.347864

8   0.759799   0.590509   0.344028   0.76854

9   0.811096   0.665767   0.217561   0.304797

10   0.0774231   0.367943   0.768124   0.0142936

11   0.941321   0.232096   0.970743   0.167941

12   0.590509   0.344028   0.05035   0.242746

13   0.665767   0.217561   0.71949   0.509156

14   0.367943   0.768124   0.739075   0.906678

15   0.232096   0.970743   0.347864   0.609261

16   0.344028   0.05035   0.76854   0.72664

17   0.217561   0.71949   0.304797   0.173989

18   0.768124   0.739075   0.0142936   0.512383

19   0.970743   0.347864   0.167941   0.0384586

 

[yasuharu@localhost CheckJumpLinux]$

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

上のプログラムはヘッダーファイルなどとともに一式を圧縮ファイルCheckJumpLinux.zipとしてまとめた。圧縮ファイル名CheckJumpLinux.zipをクリックしてダウンロードすることができる。

乱数生成のJump2toj(j)による跳び越しについては、Visual C++用のプログラム例のウェブサイトで詳しく説明した。

 

３２ビットLinux用

 

Linuxの３２ビット版では、long intは３２ビットの整数型であるが、６４ビットの整数型はlong long intになっている。このlong long int型の定数は、末尾にLではなく、LLを付けて表される。６４ビット版OS用として作成した上のプログラムを３２ビット版OS用に書き直したプログラム（３２ビット版ubuntu 12上で作成）をCheckJumpLinux32.zipにまとめた。圧縮ファイル名CheckJumpLinux32.zipをクリックしてダウンロード解凍すれば利用できる。

 

 

なお、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