next up previous
Next: 2.4 他の言語でのメルセンヌ・ツイスターの利用 Up: 2 とりあえずこれを使おう メルセンヌ・ツイスター Previous: 2.2 2002年版(mt19937ar.c)を使う


2.3 SFMT を使う

(工事中。これを学生に勧めようかと思ったけれど、ちょっと敷居が高いので、 より簡単な mt19937ar の勧め (前項 2.2) を書いた次第。)

Mersenne Twister Home Page


以下のような関数が用意されている。 

void sfmt_init_gen_rand() initializes the generator with a 32-bit integer seed.
void sfmt_init_by_array() initializes the generator with an array of 32-bit integers as the seeds.
uint32_t sfmt_genrand_uint32() generates and returns a pseudorandom 32-bit unsigned integer.
uint32_t sfmt_genrand_uint64() generates and returns a pseudorandom 64-bit unsigned integer.
void sfmt_fill_array32() fills the user-specified array with 32-bit pseudorandom integers.
void sfmt_fill_array64() fills the user-specified array with 64-bit pseudorandom integers.
inline static double sfmt_to_real1(uint32_t v) convert v to double on [0,1]-real-interval.
inline static double sfmt_genrand_real1(sfmt_t * sfmt) generates a random number on [0,1]-real-interval.
inline static double sfmt_to_real2(uint32_t v) convert v to double on [0,1)-real-interval.
inline static double sfmt_genrand_real2(sfmt_t * sfmt) generates a random number on [0,1)-real-interval.
inline static double sfmt_to_real3(uint32_t v) convert v to double on (0,1)-real-interval.
inline static double sfmt_genrand_real3(sfmt_t * sfmt) generates a random number on (0,1)-real-interval.
inline static double sfmt_to_res53(uint64_t v) convert v to double on [0,1) with 53-bit resolution
inline static double sfmt_genrand_res53(sfmt_t * sfmt) generates a random number on [0,1) with 53-bit resolution
inline static double sfmt_to_res53_mix(uint32_t x, uint32_t y) convert v to double on [0,1) with 53-bit resolution from two 32bit integers.
inline static double sfmt_genrand_res53_mix(sfmt_t * sfmt) generates a random number on [0,1) with 53-bit resolution using 32bit generation.

コンパイルの仕方とかは、 How to compile SFMT

以下は本物の gcc が使えるとした場合の話。

寄り道: gcc のインストール
MacPorts の更新をしてあれば、 例えば gcc version 6 は次のようにインストールできる (少し時間がかかる)。 
sudo port install gcc6
sudo port select --set gcc mp-gcc6

コンパイルしてチェックする
curl -O http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/SFMT-src-1.5.1.tar.gz
tar xzf SFMT-src-1.5.1.tar.gz
cd SFMT-src-1.5.1
make std>& make.log
./test-std-M19937 -b32
./test-std-M19937 -b32 > foo.txt
diff -w foo.txt SFMT.19937.out.txt

cat make.log
(diff が何も出力しなければうまく行っていると考えられる。)

open html/index.html とするとドキュメントが読める。

どうやって使うのかについては、松本氏自身のサンプル・プログラムがある (学生の好きなモンテカルロ法で円周率を求める、というやつ)。

sample1.c

/*
 * sample1.c --- a very simple program sample1.c which calculates PI using Monte-Carlo method.
 *   http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/howto-compile.html
 *   How to compile: in the directory named "SFMT-src-1.5.1"
 *      gcc -O3 -DSFMT_MEXP=607 -o sample1 SFMT.c sample1.c
 */

#include <stdio.h>
#include <stdlib.h>
#include "SFMT.h"

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

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

SFMT-src-1.5.1 というディレクトリィでコンパイル＆実行
curl -O http://nalab.mind.meiji.ac.jp/~mk/misc/20190622/sample1.c
gcc -O3 -DSFMT_MEXP=607 -o sample1 sample1.c SFMT.c
./sample1
これで 3.124800 と表示されるはず。

うまく動くようだから、では、ライブラリィを作ろう、と考えたのだが、 SFMT_MEXP の値はコンパイル時に指定する必要がある。 ということはプログラムを書き換えない限り、選択肢は次の2つ。

(a)
必要なときにその都度、 SFMT.c を -DSFMT_MEXP= を指定してコンパイルする。
(b)
例えば、デフォールトの設定 -DSFMT_MEXP=19937 でライブラリィを作り、 それをリンクして使う。

(b) の方針で行ってみよう。 最後に表示させた make.log の内容を参考にして、 次のようにライブラリィをインストールする。 
gcc -O3 -finline-functions -fomit-frame-pointer -DNDEBUG -fno-strict-aliasing \
  --param max-inline-insns-single=1800 -Wmissing-prototypes -Wall -std=c99 -c SFMT.c
ar cr libSFMT.a SFMT.o
ranlib libSFMT.a
sudo cp -p libSFMT.a /usr/local/lib
sudo mkdir /usr/local/include/SFMT
sudo cp -p *.h /usr/local/include/SFMT

インストールしたライブラリィを使ってsample1.cをコンパイル
gcc -I/usr/local/include/SFMT -O3 -DSFMT_MEXP=19937 -o sample1 sample1.c -lSFMT
./sample1
(-DSFMT_MEXP=19937 はなくてもコンパイルできるが、 警告が出る。ちょっとうっとうしい。) 

/*
 * mytest-SFMT.c
 *   gcc -I /usr/local/include/SFMT -DSFMT_MEXP=19937 mytest-SFMT.c -lSFMT
 */

#include <stdio.h>
#include <SFMT.h>

int main(void)
{
  sfmt_t sfmt;
  uint32_t seed = 1234; // 適当
  sfmt_init_gen_rand(&sfmt, seed);

  // [0,1] 上の乱数
  for (i=0; i<10; i++)
    printf("%f\n", sfmt_genrand_real1(&sfmt));

  // [0,1) 上の乱数
  for (i=0; i<10; i++)
    printf("%f\n", sfmt_genrand_real1(&sfmt));

  // [0,1) 上の乱数
  for (i=0; i<10; i++)
    printf("%f\n", sfmt_genrand_res53(&sfmt));

  return 0;
}

curl -O http://nalab.mind.meiji.ac.jp/~mk/mytest-SFMT.c


桂田 祐史