42#define MATRIX_A 0x9908b0df
43#define UPPER_MASK 0x80000000
44#define LOWER_MASK 0x7fffffff
47#define TEMPERING_MASK_B 0x9d2c5680
48#define TEMPERING_MASK_C 0xefc60000
49#define TEMPERING_SHIFT_U(y) (y >> 11)
50#define TEMPERING_SHIFT_S(y) (y << 7)
51#define TEMPERING_SHIFT_T(y) (y << 15)
52#define TEMPERING_SHIFT_L(y) (y >> 18)
54static unsigned int mt[N];
56double mersenne_array[N];
59void seed_mersenne(
unsigned long seed) {
67 mt[0] = seed & 0xffffffff;
68 for (mti = 1; mti < N; mti++)
69 mt[mti] = (69069 * mt[mti - 1]) & 0xffffffff;
79 static unsigned int mag01[2] = {0x0, MATRIX_A};
83 printf(
"DUMMY: you did not seed the generator!\n");
89 for (kk = 0; kk < N - M; kk++) {
90 y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
91 mt[kk] = mt[kk + M] ^ (y >> 1) ^ mag01[y & 0x1];
93 for (; kk < N - 1; kk++) {
94 y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
95 mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1];
97 y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
98 mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ mag01[y & 0x1];
100 for (kk = 0; kk < N; kk++) {
102 y ^= TEMPERING_SHIFT_U(y);
103 y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
104 y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
105 y ^= TEMPERING_SHIFT_L(y);
107 (double)y * 2.3283064365386963e-10;
111 return (mersenne_array[--mersenne_i]);