C言語でのメルセンヌ・ツイスタライブラリ（SFMT）の使い方

C言語の標準ライブラリに付属しているrand関数はC言語学習者であれば誰もが使用したことがあるかと思いますが、rand関数はその使い道によっては乱数としての性能が不十分な場合があったり、そもそも利用すべきでないこともあります。

本記事ではメルセンヌ・ツイスタの開発者の一人である松本眞氏により公開されているC言語向けのメルセンヌ・ツイスタライブラリ（SFMT）の利用方法を紹介します。

今すぐ使い方を知りたいという方は目次からSFMTライブラリの使い方へ飛んでください。

rand関数
- 問題１．乱数の上限
- 問題２．再現性の無さ
メルセンヌ・ツイスタ
SFMTライブラリの使い方
最後に

rand関数

C言語の標準ライブラリにあるrand関数はstdlib.hをインクルードするだけで使うことができ，非常に手軽に利用することができます。ただし、用途によってはなるべく利用しない方が良いものでもあります。

問題１．乱数の上限

まず、rand関数は環境によって2¹⁶ = 65536程度の範囲の数しか生成できないことがあります。つまり要素数7万個の配列からランダムに値を取り出すなどという場合には使えません。

さすがに最近のコンパイラでは32bitまで対応していることが多いと思いますが、64bit環境が広がっている中では32bitでも足りないという場面もありえます。

そもそも、生成される乱数の上限が環境によって異なるということ自体が問題になりえます。別のサーバーにプログラムを移したらなんだか上手く動かない…なんてことになりかねません。

問題２．再現性の無さ

また値の範囲だけでなく、生成される乱数も環境依存となります。

学術分野のように、誰にとっても実験結果が再現できるようになっていなければならない場面では、もちろん乱数も再現できなければなりません。しかしrand関数はsrand関数によって同じシード値を与えたとしても、環境によって生成される乱数が同じになるとは限りません。

また学術分野でないとしても、開発中はバグの再現のために開発者全員の環境で同じ結果が得られる方が、デバッグの効率は格段に向上します（C言語経験者であれば再現性の無いバグのツラさはよくわかるのではないでしょうか）。

メルセンヌ・ツイスタ

乱数（正確には擬似乱数）をいかに生成するかという問題は数学・アルゴリズムの分野で古くから研究されています。様々な特徴を持った乱数生成法が提案されていますが（ここでは詳しいことには触れませんが）、その中でも質の良い乱数生成方法として有名なのが メルセンヌ・ツイスタ です。松本眞氏、西村拓士氏により提案され、現在ではC++やPython、Rubyなど主要な言語でも標準で採用されている方法です。

特長１．周期が長い

擬似乱数は必ず周期性を持ち、生成するうちにいつかは最初に戻ってしまうのですが、メルセンヌ・ツイスタはその周期が2¹⁹⁹³⁷-1と非常に長く、現実的に問題となることがありません（ここで紹介するライブラリではさらに長周期にも対応しています）。

これがどれくらい長いかと言うと、例えば1ナノ秒ごとに1億台のコンピューターで同時に乱数を取得しても宇宙の寿命までに1周期が終わりません。いえ、もはやそんなレベルですらなく、100億年が10¹⁷秒、宇宙に存在する原子の数が10⁸⁰個、その一方で2³⁰⁰が10⁹⁰程度と考えると、2¹⁹⁹³⁷-1が "天文学的数字" すら余裕で飛び越えるくらいの巨大な数であることがわかります。

特長２．値間の相関が小さい

例えば連続する2つの値をそれぞれx座標とy座標の値とし、2次元平面にどんどんプロットしていく、といった使い方をする場合、連続する2つの値に相関があればプロット結果に偏りが生じます（例えば縞模様が見える、といったような）。

2次元では大丈夫としても、これが3次元、4次元…となっていくとどこかで相関が見えてきてしまうのですが、メルセンヌ・ツイスタは623次元でも均等に分布します。

つまり、上記のように座標の値として用いても問題無いですし、連続する32bit乱数を4つ繋げて128bitとして使う、といったようなことをしても乱数として十分機能します。

その他

松本眞氏によりメルセンヌ・ツイスタの特徴や名前の由来などが公開されています。論文へもアクセス可能です。

Mersenne Twister: A random number generator (since 1997/10)

またWikipediaのメルセンヌ・ツイスタのページでもわかりやすく解説されています。

メルセンヌ・ツイスタ - Wikipedia

SFMTライブラリの使い方

以下ではメルセンヌ・ツイスタの開発者の一人である松本眞氏により公開されているC言語向けのメルセンヌ・ツイスタライブラリ（SFMT）の利用方法を紹介します。

ダウンロード

まずSFMTのページの『Download SFMT which supports ...』からソースをダウンロードします。2018年2月現在、最新バージョンは1.5.1です。

テスト

次にSFMTを利用可能かどうかのチェックを行います。以下ではLinuxやMacなどUNIXコマンドが使用できる環境を前提として解説します。（詳細は How to compile SFMT を御覧ください）

ダウンロードしたファイルを解凍し、できあがったディレクトリへ移動してください。

コマンドラインから make std-check を実行します。テスト用のコードがビルドされ、出力が想定されたものであるかチェックされます。

$ make std-check
./check.sh 32 test-std
test-std-M607 output check OK
test-std-M1279 output check OK
test-std-M2281 output check OK
test-std-M4253 output check OK
test-std-M11213 output check OK
test-std-M19937 output check OK
test-std-M44497 output check OK
test-std-M86243 output check OK
test-std-M132049 output check OK
test-std-M216091 output check OK

すべて OK となっていれば大丈夫です。

次に make sse2-check を実行します。 CPUがSSE2に対応している場合に最適化してコンパイルし、その出力をチェックします。

$ make sse2-check
./check.sh 32 test-sse2
test-sse2-M607 output check OK
test-sse2-M1279 output check OK
test-sse2-M2281 output check OK
test-sse2-M4253 output check OK
test-sse2-M11213 output check OK
test-sse2-M19937 output check OK
test-sse2-M44497 output check OK
test-sse2-M86243 output check OK
test-sse2-M132049 output check OK
test-sse2-M216091 output check OK

これもすべて OK となっていればSSE2を使用できます。

ちなみに実行ファイルに -s オプションをつけて実行することで乱数生成の速度を測ることができます。

$ ./test-std-M19937 -s
32 bit BLOCK:209ms for 100000000 randoms generation
32 bit SEQUE:311ms for 100000000 randoms generation
64 bit BLOCK:223ms for 50000000 randoms generation
64 bit SEQUE:273ms for 50000000 randoms generation

$ ./test-sse2-M19937 -s
32 bit BLOCK:64ms for 100000000 randoms generation
32 bit SEQUE:143ms for 100000000 randoms generation
64 bit BLOCK:60ms for 50000000 randoms generation
64 bit SEQUE:97ms for 50000000 randoms generation

$ ./test-sse2-M607 -s
32 bit BLOCK:57ms for 100000000 randoms generation
32 bit SEQUE:129ms for 100000000 randoms generation
64 bit BLOCK:58ms for 50000000 randoms generation
64 bit SEQUE:79ms for 50000000 randoms generation

ソースコードでの利用

では実際に自分のソースコードで利用してみます。

（詳細はダウンロードしたファイルの中の html/index.html を御覧ください）

まずダウンロードしたファイルからヘッダファイル（*.h）とソースファイル（*.c）を同じディレクトリに移動させた上で、以下のように呼び出します（もちろんヘッダをインクルードしてソースを一緒にコンパイルできれば同じディレクトリである必要はありません）。

/* まずSFMT.hをインクルード */
#include "SFMT.h"

int main( int argc, char *argv[] ) {
    /* 状態を保持する構造体 */
    sfmt_t sfmt;

    /* シードを指定して初期化 */
    int seed = 0;
    sfmt_init_gen_rand( &sfmt, seed );

    /* 32bit整数を生成する場合： */
    uint32_t r_int_32 = sfmt_genrand_uint32( &sfmt );

    /* 64bit整数を生成する場合（ただし64bit環境のみ）： */
    uint64_t r_int_64 = sfmt_genrand_uint64( &sfmt );

    /* 0以上1未満の実数を生成する場合： */
    double r_real = sfmt_genrand_real2( &sfmt );

    /* 0以上1未満の52bit精度実数を生成する場合： */

    /* 1. 64bit環境ではこちら */
    uint64_t v = sfmt_genrand_uint64( &sfmt );
    double r_real_52 = sfmt_to_res53( v );

    /* 2. 32bit環境ではこちら */
    uint32_t v1 = sfmt_genrand_uint32( &sfmt );
    uint32_t v2 = sfmt_genrand_uint32( &sfmt );
    double r_real_52 = sfmt_to_res53_mix( v1, v2 );

    return 0;
}

コンパイルするときは SFMT.c をコンパイルしてリンクする必要があります。また周期を SFMT_MEXP として定義する必要があります。 SSE2を使用する場合は -msse2 を指定します。

$ gcc -O3 -msse2 -DSFMT_MEXP=19937 SFMT.c test.c

余談

標準のrand関数のように使いたい（グローバルに同じ状態の乱数でよい）場合は以下のようなコードを書けば同じように利用できます。

#include "SFMT.h"

static sfmt_t sfmt;

uint32_t my_rand() {
    return sfmt_genrand_uint32( &sfmt );
}

void my_srand( uint32_t seed ) {
    sfmt_init_gen_rand( &sfmt, seed );
}