【C】外部からヘッダーファイルをダウンロードしてMersenne Twister（MT/SFMT）を利用する
Cでもメルセンヌツイスターを利用できます。そんなメモ。
はじめに
以前、C++11でMersenne Twisterを利用する記事を書きました。
dandy-tech.hatenablog.jp
すると今度は「C言語ではMersenne Twisterを利用するにはどうすれば...。」という疑問が出てきます。
結論を言うと外部から必要ファイルのダウンロードをするだけでOK。
まあ今時そんなもんですよね。
ファイルのダウンロード
様々なサイトでファイルがダウンロードできるようです。
www.math.sci.hiroshima-u.ac.jp
速度が必要な場合にはMTではなく以下のSFMTの方がずっと好ましいです。
SIMD-oriented Fast Mersenne Twister (SFMT)
www.math.sci.hiroshima-u.ac.jp
速度の比較についてはこちら
せっかくなので、MTとSFMTの双方の利用法を記述します。
僕の経験ではMTの方がコンパイルに余計なオプションが必要ないので、手軽に利用できる印象です。
MTの場合
今回は64bit版をダウンロードしました。
$ wget http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt19937-64.tgz $ tar zxvf tar zxvf mt19937-64.tgz
SFMTの場合
記事執筆時はSFMT ver1.5.1が最新でしたので、こちらをダウンロード＆解凍。
$ wget http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/SFMT-src-1.5.1.tar.gz $ tar zxvf SFMT-src-1.5.1.tar.gz
コンパイル＆実行
MTの場合
mt64.hの中に乱数の種類ごとの関数が記述されています。
例えば
// 0以上、2^64-1以下の整数乱数 /* generates a random number on [0, 2^64-1]-interval */ unsigned long long genrand64_int64(void); // 0以上、1未満の実数乱数 /* generates a random number on [0,1)-real-interval */ double genrand64_real2(void);
などです。
試しに以下のような条件で一次元ランダムウォークさせてみます。
「1000回サイコロを振る。サイコロの出目（１～６）だけx軸方向を移動する。ただし、サイコロを振った回数が奇数回目なら+x方向、偶数回目なら-x方向に移動すること。」
#include <stdio.h> #include <stdlib.h> #include <time.h> #include "mt19937-64/mt64.h" int main(void) { int i, dice, x; int p = 0; //position int N = 1000 time_t t; t = time(0); init_genrand64(t); //mtをinitialize FILE *fp; fp = fopen("random_walk.dat","w"); for (i=1; i<N+1; i++) { dice = genrand64_int64() % 6 + 1; if(i%2==0){ x = -1 * dice; // printf("%2d\n", x); }else{ x = dice; // printf("%2u \n", dice); } p = p + x; if(i==N){ fprintf(fp, "%4d %2d", i, p); break; } fprintf(fp, "%4d %2d\n", i, p); } fclose(fp); return 0; }
ここではtime.hをincludeして乱数のseedに用いています。
$ gcc random_walk.c mt19937-64/mt19937-64.c $ ./a.out
SFMTの場合
ダウンロードページの最下部までスクロールするとコンパイル方法のページのリンクがあります。
念のためリンクを貼っておきます。
（解凍したファイルの中にもhtmlファイルがありますが...）
How to compile SFMT
このHow toでは解凍したファイル内のtest.cを例にコンパイル＆実行がなされています。
SFMTについても利用できる関数はSFMT.hを参照するといいでしょう。
同様の条件のランダムウォークを生成するファイルはこちら
#include <stdio.h> #include <stdlib.h> #include <time.h> #include "SFMT-src-1.5.1/SFMT.h" int main(void) { int i, dice, x; int p = 0; //position int N = 1000; time_t t; sfmt_t sfmt; t = time(0); sfmt_init_gen_rand(&sfmt, t); // printf("%"PRIu64 "\n", sfmt_genrand_uint64(&sfmt)); FILE *fp; fp = fopen("random_walk.dat","w"); for (i=1; i<N+1; i++) { dice = sfmt_genrand_uint64(&sfmt) % 6 + 1; if(i%2==0){ x = -1 * dice; // printf("%2d\n", x); }else{ x = dice; // printf("%2u \n", dice); } p = p + x; printf("%2d\n", p); if(i==N){ fprintf(fp, "%4d %2d", i, p); break; } fprintf(fp, "%4d %2d\n", i, p); } fclose(fp); return 0; }
$ gcc -O3 -DSFMT_MEXP=19937 SFMT-src-1.5.1/SFMT.c random_walk2.c $ ./a.out
描画してみる
今回もROOT６で横軸をサイコロを振った回数N、縦軸をxとしてプロットしました。
