2.2 2002年版(mt19937ar.c)を使う

2019年現在、松本氏は改良版である SFMT を使うことを推奨しているようだが、 オリジナルに近い 2002年版は使いやすいと思う (この項と次項 2.3 を見比べてみよう)。

この項は、以前書いた 「疑似乱数についてのメモ」 (http://nalab.mind.meiji.ac.jp/~mk/labo/text/giji-ransuu/) の付録部分を加筆修正したものである。

本家 「Mersenne Twister Home Page」 の C プログラム mt19937ar.sep.tgz を試してみたので紹介する。

必要な関数の入っている mt19937ar.c をコンパイルして、 オブジェクト・ファイルを作り、テスト・プログラム mtTest.c をコンパイル＆リンクし、 実行結果を付属する mt19937ar.out と比較する。

最初の行以外は、ターミナルにそのままコピペして動くはず

cd どこか適当なディレクトリィ
curl -O http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.sep.tgz
tar xzf mt19937ar.sep.tgz
gcc -O -c mt19937ar.c
gcc -o mtTest mtTest.c mt19937ar.o
./mtTest > foo
diff foo mt19937ar.out

これで何も表示されなければ成功ということ。 私は本物の gcc を使ったが、macOS Xcode の cc や gcc でも大丈夫みたいだ。

インストールするには、少々強引だが、

続いてこうやってインストール

ar cr libmt19937ar.a mt19937ar.o
ranlib libmt19937ar.a
sudo cp -p libmt19937ar.a /usr/local/lib
sudo cp -p mt19937ar.h /usr/local/include

(sudo でパスワードを尋ねられるので、 最後の1行は別にコピペするべきかも。)

これで C プログラムの中に

#include <mt19937ar.h>

として、 リンクするときに -lmt19937ar というリンカー・オプションを指定する。

次に掲げるのは、有名なモンテカルロ法で円周率を求める、というプログラムである。

sample1-mt19937ar.c

/*
 * sample1.c --- a very simple program which calculates PI using Monte-Carlo method.
 *   http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/howto-compile.html
 *   にある sample1.c (SFMT用) を mt19937ar 用に書き換えたもの
 *      gcc -O3 -o sample1 sample1.c -lmt19937ar
 */

#include <stdio.h>
#include <stdlib.h>
#include <mt19937ar.h>

int main(int argc, char *argv[])
{
    int i, cnt;
    unsigned long seed;
    double x, y, pi;
    const int NUM = 10000;

    if (argc >= 2) {
        seed = strtol(argv[1], NULL, 10);
    } else {
        seed = 12345;
    }
    cnt = 0;
    init_genrand(seed);
    for (i = 0; i < NUM; i++) {
        x = genrand_real2();
        y = genrand_real2();
        if (x * x + y * y < 1.0) {
            cnt++;
        }
    }
    pi = (double) cnt / NUM * 4;
    printf("%lf\n", pi);
    return 0;
}

こうやって入手・コンパイル・実行

curl -O http://nalab.mind.meiji.ac.jp/~mk/misc/20190622/sample1-mt19937ar.c
gcc -I/usr/local/include sample1-mt19937ar.c -lmt19937ar
./a.out

これで 3.132400 と表示されるはず。

初期化は

  unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
  init_by_array(init, length);

のように unsigned long int の配列データで初期化するか、

  init_genrand((unsigned long)20120328);

のように unsigned long int 型データで初期化する。 もちろん (srand() でやったように)

  init_genrand((unsigned long)time(NULL));

のように time() を使っても良いだろう。

• unsigned long genrand_int32(void) 0 から 0xffffffff まで
• long genrand_int31(void) 0 から $2^{31}-1$ まで。
• double genrand_real1(void) $[0,1]$
• double genrand_real2(void) $[0,1)$
• double genrand_real3(void) $(0,1)$
• double genrand_res53(void) $[0,1)$

C++ では、新しい規格ではメルセンヌ・ツイスターが含まれているそうだが (もう g++ で使えるの？)、 この C バージョンも利用可能である。

(1)

C++ プログラム中に

extern "C" {
#include "mt19937ar.h"
};

のようにインクルードする。

(2)

mt19937ar.c を C++ プログラムとして利用する。

cp -p mt199e7ar.c mt19937ar.cpp
g++ -O -c mt19937ar.cpp