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
|
#include "mt.h"
#include <stdlib.h>
#include <stdio.h>
/* This code is based on mt19937ar.c, written by Takuji Nishimura and
Makoto Matsumoto (20020126). Further details are available at
<http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html>.
REFERENCE
M. Matsumoto and T. Nishimura,
"Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
Pseudo-Random Number Generator",
ACM Transactions on Modeling and Computer Simulation,
Vol. 8, No. 1, January 1998, pp 3--30. */
void mt_init_seed( struct mt *m, uint32_t seed )
{
int i;
uint32_t *mt;
mt = m->mt;
mt[0] = seed & 0xffffffff;
for ( i = 1; i < N; i++ )
mt[i] = 1812433253 * (mt[i-1]^(mt[i-1]>>30)) + i;
m->mti = N;
}
struct mt *mt_setup(uint32_t seed)
{
struct mt *self = malloc(sizeof(struct mt));
if (self)
mt_init_seed( self, seed );
return self;
}
struct mt *mt_setup_array( uint32_t *array, int n )
{
int i, j, k;
struct mt *self = malloc(sizeof(struct mt));
uint32_t *mt;
if (self) {
mt_init_seed( self, 19650218UL );
i = 1; j = 0;
k = ( N > n ? N : n );
mt = self->mt;
for (; k; k--) {
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
+ array[j] + j;
i++; j++;
if (i>=N) { mt[0] = mt[N-1]; i=1; }
if (j>=n) j=0;
}
for (k=N-1; k; k--) {
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) - i;
i++;
if (i>=N) { mt[0] = mt[N-1]; i=1; }
}
mt[0] = 0x80000000UL;
}
return self;
}
void mt_free(struct mt *self)
{
free(self);
}
/* Returns a pseudorandom number which is uniformly distributed in [0,1) */
double mt_genrand(struct mt *self)
{
int kk;
uint32_t y;
static uint32_t mag01[2] = {0x0, 0x9908b0df};
static const uint32_t UP_MASK = 0x80000000, LOW_MASK = 0x7fffffff;
if (self->mti >= N) {
for (kk = 0; kk < N-M; kk++) {
y = (self->mt[kk] & UP_MASK) | (self->mt[kk+1] & LOW_MASK);
self->mt[kk] = self->mt[kk+M] ^ (y >> 1) ^ mag01[y & 1];
}
for (; kk < N-1; kk++) {
y = (self->mt[kk] & UP_MASK) | (self->mt[kk+1] & LOW_MASK);
self->mt[kk] = self->mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 1];
}
y = (self->mt[N-1] & UP_MASK) | (self->mt[0] & LOW_MASK);
self->mt[N-1] = self->mt[M-1] ^ (y >> 1) ^ mag01[y & 1];
self->mti = 0;
}
y = self->mt[self->mti++];
y ^= y >> 11;
y ^= y << 7 & 0x9d2c5680;
y ^= y << 15 & 0xefc60000;
y ^= y >> 18;
return y*(1.0/4294967296.0);
}
|
