1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
|
// This file is the C Code from the Mersenne Twister wikipedia article
// https://en.wikipedia.org/wiki/Mersenne_Twister#C_code
// bodged together with the init_by_array and main functions from
// https://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/MT2002/CODES/mt19937ar.c
// so I can test if the C Code on wikipedia generates the correct random numbers
// as seen in
// https://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/MT2002/CODES/mt19937ar.out
//
// I did not check them all, but a few at random. They match, which makes me think
// they operate the same. Now I need to look for the bug in my Rust,,,
#include <stdint.h>
#include <stdio.h>
#define n 624
#define m 397
#define w 32
#define r 31
#define UMASK (0xffffffffUL << r)
#define LMASK (0xffffffffUL >> (w-r))
#define a 0x9908b0dfUL
#define u 11
#define s 7
#define t 15
#define l 18
#define b 0x9d2c5680UL
#define c 0xefc60000UL
#define f 1812433253UL
typedef struct
{
uint32_t state_array[n]; // the array for the state vector
int state_index; // index into state vector array, 0 <= state_index <= n-1 always
} mt_state;
void initialize_state(mt_state* state, uint32_t seed)
{
uint32_t* state_array = &(state->state_array[0]);
state_array[0] = seed; // suggested initial seed = 19650218UL
for (int i=1; i<n; i++)
{
seed = f * (seed ^ (seed >> (w-2))) + i; // Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
state_array[i] = seed;
}
state->state_index = 0;
}
/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
/* slight change for C++, 2004/2/26 */
void init_by_array(mt_state* state, unsigned long init_key[], int key_length)
{
int i, j, k;
initialize_state(state, 19650218UL);
i=1; j=0;
k = (n>key_length ? n : key_length);
uint32_t* mt = &(state->state_array[0]);
for (; k; k--) {
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
+ init_key[j] + j; /* non linear */
mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
i++; j++;
if (i>=n) { mt[0] = mt[n-1]; i=1; }
if (j>=key_length) j=0;
}
for (k=n-1; k; k--) {
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
- i; /* non linear */
mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
i++;
if (i>=n) { mt[0] = mt[n-1]; i=1; }
}
mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
}
uint32_t random_uint32(mt_state* state)
{
uint32_t* state_array = &(state->state_array[0]);
int k = state->state_index; // point to current state location
// 0 <= state_index <= n-1 always
// int k = k - n; // point to state n iterations before
// if (k < 0) k += n; // modulo n circular indexing
// the previous 2 lines actually do nothing
// for illustration only
int j = k - (n-1); // point to state n-1 iterations before
if (j < 0) j += n; // modulo n circular indexing
uint32_t x = (state_array[k] & UMASK) | (state_array[j] & LMASK);
uint32_t xA = x >> 1;
if (x & 0x00000001UL) xA ^= a;
j = k - (n-m); // point to state n-m iterations before
if (j < 0) j += n; // modulo n circular indexing
x = state_array[j] ^ xA; // compute next value in the state
state_array[k++] = x; // update new state value
if (k >= n) k = 0; // modulo n circular indexing
state->state_index = k;
uint32_t y = x ^ (x >> u); // tempering
y = y ^ ((y << s) & b);
y = y ^ ((y << t) & c);
uint32_t z = y ^ (y >> l);
return z;
}
int main(void)
{
int i;
unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
mt_state state;
state.state_index = 0;
initialize_state(&state, 19650218UL);
//init_by_array(&state, init, length);
/*printf("1000 outputs of genrand_int32()\n");
for (i=0; i<1000; i++) {
printf("%10lu ", random_uint32(&state));
if (i%5==4) printf("\n");
}*/
// Print the internal state so we can check it against the Rust
for (i = 0; i < n; ++i) {
printf("%10lu\n", state.state_array[i]);
}
return 0;
}
|
