44#define MATRIX_A 0x9908b0df
45#define UPPER_MASK 0x80000000
46#define LOWER_MASK 0x7fffffff
49#define TEMPERING_MASK_B 0x9d2c5680
50#define TEMPERING_MASK_C 0xefc60000
51#define TEMPERING_SHIFT_U(y) (y >> 11)
52#define TEMPERING_SHIFT_S(y) (y << 7)
53#define TEMPERING_SHIFT_T(y) (y << 15)
54#define TEMPERING_SHIFT_L(y) (y >> 18)
56static unsigned int mt[N];
57static int mersenne_i = -1;
60void seed_mersenne(
long seed) {
68 mt[0] = seed & 0xffffffff;
69 for (mti = 1; mti < N; mti++)
70 mt[mti] = (69069 * mt[mti - 1]) & 0xffffffff;
80 static unsigned int mag01[2] = {0x0, MATRIX_A};
84 printf(
"MERSENNE: you did not seed the generator!\n");
90 for (kk = 0; kk < N - M; kk++) {
91 y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
92 mt[kk] = mt[kk + M] ^ (y >> 1) ^ mag01[y & 0x1];
94 for (; kk < N - 1; kk++) {
95 y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
96 mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1];
98 y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
99 mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ mag01[y & 0x1];
103 if (mersenne_i >= N) {
107 unsigned int y = mt[mersenne_i++];
108 y ^= TEMPERING_SHIFT_U(y);
109 y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
110 y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
111 y ^= TEMPERING_SHIFT_L(y);
112 return (
double)y * 2.3283064365386963e-10;