mt.c

#include <stdio.h>

/* Period parameters */
#define N 624
#define M 397
#define MATRIX_A 0x9908b0dfUL   /* constant vector a */
#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
#define LOWER_MASK 0x7fffffffUL /* least significant r bits */

static unsigned long mt[N]; /* the array for the state vector  */
static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */

void init_genrand(unsigned long s)
{
    mt[0]= s & 0xffffffffUL;
    for (mti=1; mti<N; mti++) {
        mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); 
        mt[mti] &= 0xffffffffUL;
    }
}

unsigned long temper(unsigned long y)
{
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);
    return y;
}

/*
genrand_int32の実装3通り
それぞれで変数mtiの意味が違ってくるので注意
*/

/* ナイーブな実装 (数列の最後のN項を順番にmtに格納しておく方法) */
unsigned long genrand_int32_naive(void)
{
    unsigned long y;
    static unsigned long mag01[2]={0x0UL, MATRIX_A};
    unsigned long last = mt[0];
    /* 配列の要素を一つ左へずらす */
    for (int i = 0; i < N - 1; i ++) {
        mt[i] = mt[i+1];
    }
    y = (last&UPPER_MASK)|(mt[0]&LOWER_MASK);
    mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];

    return temper(mt[N-1]);
}

/* ラウンドロビン法 （ずらす代わりにmt[mti-k]に数列の最後から数えてk番目の値をいれる） */
unsigned long genrand_int32_roundrobin(void)
{
    unsigned long y;
    static unsigned long mag01[2]={0x0UL, MATRIX_A};
    y = (mt[mti % N]&UPPER_MASK)|(mt[(mti+1)%N]&LOWER_MASK);
    y = mt[mti % N] = mt[(mti+M)%N] ^ (y >> 1) ^ mag01[y & 0x1UL];
    mti = (mti + 1) % N;
    return temper(y);
}

/* オリジナルの実装 (N回に一回テーブルを更新する) */
unsigned long genrand_int32_orig(void)
{
    unsigned long y;
    static unsigned long mag01[2]={0x0UL, MATRIX_A};

    if (mti >= N) {
        int kk;

        for (kk=0;kk<N-M;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
        }
        for (;kk<N-1;kk++) {
            y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
            mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
        }
        y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
        mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];

        mti = 0;
    }
  
    y = mt[mti++];
    return temper(y);
}

/* どの実装でも同じ結果になることを確認 */
int main(void)
{
    int n = 10;
    init_genrand(0);
    for (int i=0; i<n; i++) {
      printf("%08x ", genrand_int32_naive());
    }
    printf("\n");
    init_genrand(0);
    for (int i=0; i<n; i++) {
      printf("%08x ", genrand_int32_roundrobin());
    }
    printf("\n");
    init_genrand(0);
    for (int i=0; i<n; i++) {
      printf("%08x ", genrand_int32_orig());
    }
    printf("\n");
    return 0;
}