2 /* This code taken from the GNU Scientific Library and adapted for use in
3 * fastq-tools. It is subject to the following licence.
6 /* This program is free software; you can redistribute it and/or
7 modify it under the terms of the GNU General Public License as
8 published by the Free Software Foundation; either version 3 of the
9 License, or (at your option) any later version.
11 This program is distributed in the hope that it will be useful, but
12 WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 General Public License for more details. You should have received
15 a copy of the GNU General Public License along with this program;
16 if not, write to the Free Foundation, Inc., 59 Temple Place, Suite
17 330, Boston, MA 02111-1307 USA
19 Original implementation was copyright (C) 1997 Makoto Matsumoto and
20 Takuji Nishimura. Coded by Takuji Nishimura, considering the
21 suggestions by Topher Cooper and Marc Rieffel in July-Aug. 1997, "A
22 C-program for MT19937: Integer version (1998/4/6)"
24 This implementation copyright (C) 1998 Brian Gough. I reorganized
25 the code to use the module framework of GSL. The license on this
26 implementation was changed from LGPL to GPL, following paragraph 3
27 of the LGPL, version 2.
31 The seeding procedure has been updated to match the 10/99 release
36 The seeding procedure has been updated again to match the 2002
39 The original code included the comment: "When you use this, send an
40 email to: matumoto@math.keio.ac.jp with an appropriate reference to
43 Makoto Matsumoto has a web page with more information about the
44 generator, http://www.math.keio.ac.jp/~matumoto/emt.html.
46 The paper below has details of the algorithm.
48 From: Makoto Matsumoto and Takuji Nishimura, "Mersenne Twister: A
49 623-dimensionally equidistributerd uniform pseudorandom number
50 generator". ACM Transactions on Modeling and Computer Simulation,
51 Vol. 8, No. 1 (Jan. 1998), Pages 3-30
53 You can obtain the paper directly from Makoto Matsumoto's web page.
55 The period of this generator is 2^{19937} - 1.
66 unsigned long int default_seed = 4357;
69 #define N 624 /* Period parameters */
72 /* most significant w-r bits */
73 static const unsigned long UPPER_MASK = 0x80000000UL;
75 /* least significant r bits */
76 static const unsigned long LOWER_MASK = 0x7fffffffUL;
79 static const unsigned long RNG_MIN = 0;
80 static const unsigned long RNG_MAX = 0xffffffffUL;
90 static inline unsigned long mt_get(rng_t* state)
93 unsigned long int *const mt = state->mt;
95 #define MAGIC(y) (((y)&0x1) ? 0x9908b0dfUL : 0)
98 { /* generate N words at one time */
101 for (kk = 0; kk < N - M; kk++)
103 unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
104 mt[kk] = mt[kk + M] ^ (y >> 1) ^ MAGIC(y);
106 for (; kk < N - 1; kk++)
108 unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);
109 mt[kk] = mt[kk + (M - N)] ^ (y >> 1) ^ MAGIC(y);
113 unsigned long y = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
114 mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ MAGIC(y);
124 k ^= (k << 7) & 0x9d2c5680UL;
125 k ^= (k << 15) & 0xefc60000UL;
133 double mt_get_double(rng_t* state)
135 return mt_get(state) / 4294967296.0;
138 void mt_set(rng_t* state, unsigned long int s)
143 s = 4357; /* the default seed is 4357 */
145 state->mt[0]= s & 0xffffffffUL;
147 for (i = 1; i < N; i++)
149 /* See Knuth's "Art of Computer Programming" Vol. 2, 3rd
150 Ed. p.106 for multiplier. */
153 (1812433253UL * (state->mt[i-1] ^ (state->mt[i-1] >> 30)) + i);
155 state->mt[i] &= 0xffffffffUL;
161 #if 0 // these two seeding procedures are not used
162 static void mt_1999_set(rng_t* state, unsigned long int s)
167 s = 4357; /* the default seed is 4357 */
169 /* This is the October 1999 version of the seeding procedure. It
170 was updated by the original developers to avoid the periodicity
171 in the simple congruence originally used.
173 Note that an ANSI-C unsigned long integer arithmetic is
174 automatically modulo 2^32 (or a higher power of two), so we can
175 safely ignore overflow. */
177 #define LCG(x) ((69069 * x) + 1) &0xffffffffUL
179 for (i = 0; i < N; i++)
181 state->mt[i] = s & 0xffff0000UL;
183 state->mt[i] |= (s &0xffff0000UL) >> 16;
190 /* This is the original version of the seeding procedure, no longer
191 used but available for compatibility with the original MT19937. */
193 static void mt_1998_set(rng_t* state, unsigned long int s)
198 s = 4357; /* the default seed is 4357 */
200 state->mt[0] = s & 0xffffffffUL;
202 #define LCG1998(n) ((69069 * n) & 0xffffffffUL)
204 for (i = 1; i < N; i++)
205 state->mt[i] = LCG1998 (state->mt[i - 1]);
211 rng_t* fastq_rng_alloc()
213 rng_t* rng = malloc_or_die(sizeof(rng_t));
214 mt_set(rng, default_seed);
219 void fastq_rng_free(rng_t* rng)
224 void fastq_rng_seed(rng_t* rng, unsigned long seed)
229 unsigned long fastq_rng_uniform_int(rng_t* rng, unsigned long k)
231 unsigned long scale = (RNG_MAX - RNG_MIN) / k;
236 r = (mt_get(rng) - RNG_MIN) / scale;