1 /* This program is free software; you can redistribute it and/or
2 modify it under the terms of the GNU General Public License as
3 published by the Free Software Foundation; either version 3 of the
4 License, or (at your option) any later version.
6 This program is distributed in the hope that it will be useful, but
7 WITHOUT ANY WARRANTY; without even the implied warranty of
8 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
9 General Public License for more details. You should have received
10 a copy of the GNU General Public License along with this program;
11 if not, write to the Free Foundation, Inc., 59 Temple Place, Suite
12 330, Boston, MA 02111-1307 USA
14 Original implementation was copyright (C) 1997 Makoto Matsumoto and
15 Takuji Nishimura. Coded by Takuji Nishimura, considering the
16 suggestions by Topher Cooper and Marc Rieffel in July-Aug. 1997, "A
17 C-program for MT19937: Integer version (1998/4/6)"
19 This implementation copyright (C) 1998 Brian Gough. I reorganized
20 the code to use the module framework of GSL. The license on this
21 implementation was changed from LGPL to GPL, following paragraph 3
22 of the LGPL, version 2.
26 The seeding procedure has been updated to match the 10/99 release
31 The seeding procedure has been updated again to match the 2002
34 The original code included the comment: "When you use this, send an
35 email to: matumoto@math.keio.ac.jp with an appropriate reference to
38 Makoto Matsumoto has a web page with more information about the
39 generator, http://www.math.keio.ac.jp/~matumoto/emt.html.
41 The paper below has details of the algorithm.
43 From: Makoto Matsumoto and Takuji Nishimura, "Mersenne Twister: A
44 623-dimensionally equidistributerd uniform pseudorandom number
45 generator". ACM Transactions on Modeling and Computer Simulation,
46 Vol. 8, No. 1 (Jan. 1998), Pages 3-30
48 You can obtain the paper directly from Makoto Matsumoto's web page.
50 The period of this generator is 2^{19937} - 1.
56 #include <gsl/gsl_rng.h>
58 static inline unsigned long int mt_get (void *vstate);
59 static double mt_get_double (void *vstate);
60 static void mt_set (void *state, unsigned long int s);
62 #define N 624 /* Period parameters */
65 /* most significant w-r bits */
66 static const unsigned long UPPER_MASK = 0x80000000UL;
68 /* least significant r bits */
69 static const unsigned long LOWER_MASK = 0x7fffffffUL;
78 static inline unsigned long
81 mt_state_t *state = (mt_state_t *) vstate;
84 unsigned long int *const mt = state->mt;
86 #define MAGIC(y) (((y)&0x1) ? 0x9908b0dfUL : 0)
89 { /* generate N words at one time */
92 for (kk = 0; kk < N - M; kk++)
94 unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
95 mt[kk] = mt[kk + M] ^ (y >> 1) ^ MAGIC(y);
97 for (; kk < N - 1; kk++)
99 unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
100 mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ MAGIC(y);
104 unsigned long y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
105 mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ MAGIC(y);
115 k ^= (k << 7) & 0x9d2c5680UL;
116 k ^= (k << 15) & 0xefc60000UL;
125 mt_get_double (void * vstate)
127 return mt_get (vstate) / 4294967296.0 ;
131 mt_set (void *vstate, unsigned long int s)
133 mt_state_t *state = (mt_state_t *) vstate;
137 s = 4357; /* the default seed is 4357 */
139 state->mt[0]= s & 0xffffffffUL;
141 for (i = 1; i < N; i++)
143 /* See Knuth's "Art of Computer Programming" Vol. 2, 3rd
144 Ed. p.106 for multiplier. */
147 (1812433253UL * (state->mt[i-1] ^ (state->mt[i-1] >> 30)) + i);
149 state->mt[i] &= 0xffffffffUL;
156 mt_1999_set (void *vstate, unsigned long int s)
158 mt_state_t *state = (mt_state_t *) vstate;
162 s = 4357; /* the default seed is 4357 */
164 /* This is the October 1999 version of the seeding procedure. It
165 was updated by the original developers to avoid the periodicity
166 in the simple congruence originally used.
168 Note that an ANSI-C unsigned long integer arithmetic is
169 automatically modulo 2^32 (or a higher power of two), so we can
170 safely ignore overflow. */
172 #define LCG(x) ((69069 * x) + 1) &0xffffffffUL
174 for (i = 0; i < N; i++)
176 state->mt[i] = s & 0xffff0000UL;
178 state->mt[i] |= (s &0xffff0000UL) >> 16;
185 /* This is the original version of the seeding procedure, no longer
186 used but available for compatibility with the original MT19937. */
189 mt_1998_set (void *vstate, unsigned long int s)
191 mt_state_t *state = (mt_state_t *) vstate;
195 s = 4357; /* the default seed is 4357 */
197 state->mt[0] = s & 0xffffffffUL;
199 #define LCG1998(n) ((69069 * n) & 0xffffffffUL)
201 for (i = 1; i < N; i++)
202 state->mt[i] = LCG1998 (state->mt[i - 1]);
207 static const gsl_rng_type mt_type =
208 {"mt19937", /* name */
209 0xffffffffUL, /* RAND_MAX */
216 static const gsl_rng_type mt_1999_type =
217 {"mt19937_1999", /* name */
218 0xffffffffUL, /* RAND_MAX */
225 static const gsl_rng_type mt_1998_type =
226 {"mt19937_1998", /* name */
227 0xffffffffUL, /* RAND_MAX */
234 const gsl_rng_type *gsl_rng_mt19937 = &mt_type;
235 const gsl_rng_type *gsl_rng_mt19937_1999 = &mt_1999_type;
236 const gsl_rng_type *gsl_rng_mt19937_1998 = &mt_1998_type;
238 /* MT19937 is the default generator, so define that here too */
240 const gsl_rng_type *gsl_rng_default = &mt_type;
241 unsigned long int gsl_rng_default_seed = 0;