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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
|
/* rng/taus.c
*
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or (at
* your option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
/* This is a maximally equidistributed combined Tausworthe
generator. The sequence is,
x_n = (s1_n ^ s2_n ^ s3_n)
s1_{n+1} = (((s1_n & 4294967294) <<12) ^ (((s1_n <<13) ^ s1_n) >>19))
s2_{n+1} = (((s2_n & 4294967288) << 4) ^ (((s2_n << 2) ^ s2_n) >>25))
s3_{n+1} = (((s3_n & 4294967280) <<17) ^ (((s3_n << 3) ^ s3_n) >>11))
computed modulo 2^32. In the three formulas above '^' means
exclusive-or (C-notation), not exponentiation. Note that the
algorithm relies on the properties of 32-bit unsigned integers (it
is formally defined on bit-vectors of length 32). I have added a
bitmask to make it work on 64 bit machines.
We initialize the generator with s1_1 .. s3_1 = s_n MOD m, where
s_n = (69069 * s_{n-1}) mod 2^32, and s_0 = s is the user-supplied
seed.
The theoretical value of x_{10007} is 2733957125. The subscript
10007 means (1) seed the generator with s=1 (2) do six warm-up
iterations, (3) then do 10000 actual iterations.
The period of this generator is about 2^88.
From: P. L'Ecuyer, "Maximally Equidistributed Combined Tausworthe
Generators", Mathematics of Computation, 65, 213 (1996), 203--213.
This is available on the net from L'Ecuyer's home page,
http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme.ps
ftp://ftp.iro.umontreal.ca/pub/simulation/lecuyer/papers/tausme.ps
Update: April 2002
There is an erratum in the paper "Tables of Maximally
Equidistributed Combined LFSR Generators", Mathematics of
Computation, 68, 225 (1999), 261--269:
http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme2.ps
... the k_j most significant bits of z_j must be non-
zero, for each j. (Note: this restriction also applies to the
computer code given in [4], but was mistakenly not mentioned in
that paper.)
This affects the seeding procedure by imposing the requirement
s1 > 1, s2 > 7, s3 > 15.
The generator taus2 has been added to satisfy this requirement.
The original taus generator is unchanged.
Update: November 2002
There was a bug in the correction to the seeding procedure for s2.
It affected the following seeds 254679140 1264751179 1519430319
2274823218 2529502358 3284895257 3539574397 (s2 < 8).
*/
static inline unsigned long int taus_get (void *vstate);
static double taus_get_double (void *vstate);
static void taus_set (void *state, unsigned long int s);
typedef struct
{
unsigned long int s1, s2, s3;
}
taus_state_t;
static inline unsigned long
taus_get (void *vstate)
{
taus_state_t *state = (taus_state_t *) vstate;
#define MASK 0xffffffffUL
#define TAUSWORTHE(s,a,b,c,d) (((s &c) <<d) &MASK) ^ ((((s <<a) &MASK)^s) >>b)
state->s1 = TAUSWORTHE (state->s1, 13, 19, 4294967294UL, 12);
state->s2 = TAUSWORTHE (state->s2, 2, 25, 4294967288UL, 4);
state->s3 = TAUSWORTHE (state->s3, 3, 11, 4294967280UL, 17);
return (state->s1 ^ state->s2 ^ state->s3);
}
static double
taus_get_double (void *vstate)
{
return taus_get (vstate) / 4294967296.0 ;
}
static void
taus_set (void *vstate, unsigned long int s)
{
taus_state_t *state = (taus_state_t *) vstate;
if (s == 0)
s = 1; /* default seed is 1 */
#define LCG(n) ((69069 * n) & 0xffffffffUL)
state->s1 = LCG (s);
state->s2 = LCG (state->s1);
state->s3 = LCG (state->s2);
/* "warm it up" */
taus_get (state);
taus_get (state);
taus_get (state);
taus_get (state);
taus_get (state);
taus_get (state);
return;
}
static void
taus2_set (void *vstate, unsigned long int s)
{
taus_state_t *state = (taus_state_t *) vstate;
if (s == 0)
s = 1; /* default seed is 1 */
#define LCG(n) ((69069 * n) & 0xffffffffUL)
state->s1 = LCG (s);
if (state->s1 < 2) state->s1 += 2UL;
state->s2 = LCG (state->s1);
if (state->s2 < 8) state->s2 += 8UL;
state->s3 = LCG (state->s2);
if (state->s3 < 16) state->s3 += 16UL;
/* "warm it up" */
taus_get (state);
taus_get (state);
taus_get (state);
taus_get (state);
taus_get (state);
taus_get (state);
return;
}
static const gsl_rng_type taus_type =
{"taus", /* name */
0xffffffffUL, /* RAND_MAX */
0, /* RAND_MIN */
sizeof (taus_state_t),
&taus_set,
&taus_get,
&taus_get_double};
const gsl_rng_type *gsl_rng_taus = &taus_type;
static const gsl_rng_type taus2_type =
{"taus2", /* name */
0xffffffffUL, /* RAND_MAX */
0, /* RAND_MIN */
sizeof (taus_state_t),
&taus2_set,
&taus_get,
&taus_get_double};
const gsl_rng_type *gsl_rng_taus2 = &taus2_type;
|
