| /* rng/tt.c |
| * |
| * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 2 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 the TT800 twisted GSFR generator for 32 bit integers. It |
| has been superceded by MT19937 (mt.c). The period is 2^800. |
| |
| This implementation is based on tt800.c, July 8th 1996 version by |
| M. Matsumoto, email: matumoto@math.keio.ac.jp |
| |
| From: Makoto Matsumoto and Yoshiharu Kurita, "Twisted GFSR |
| Generators II", ACM Transactions on Modelling and Computer |
| Simulation, Vol. 4, No. 3, 1994, pages 254-266. */ |
| |
| static inline unsigned long int tt_get (void *vstate); |
| static double tt_get_double (void *vstate); |
| static void tt_set (void *state, unsigned long int s); |
| |
| #define N 25 |
| #define M 7 |
| |
| typedef struct |
| { |
| int n; |
| unsigned long int x[N]; |
| } |
| tt_state_t; |
| |
| static inline unsigned long int |
| tt_get (void *vstate) |
| { |
| tt_state_t *state = (tt_state_t *) vstate; |
| |
| /* this is the magic vector, a */ |
| |
| const unsigned long mag01[2] = |
| {0x00000000, 0x8ebfd028UL}; |
| unsigned long int y; |
| unsigned long int *const x = state->x; |
| int n = state->n; |
| |
| if (n >= N) |
| { |
| int i; |
| for (i = 0; i < N - M; i++) |
| { |
| x[i] = x[i + M] ^ (x[i] >> 1) ^ mag01[x[i] % 2]; |
| } |
| for (; i < N; i++) |
| { |
| x[i] = x[i + (M - N)] ^ (x[i] >> 1) ^ mag01[x[i] % 2]; |
| }; |
| n = 0; |
| } |
| |
| y = x[n]; |
| y ^= (y << 7) & 0x2b5b2500UL; /* s and b, magic vectors */ |
| y ^= (y << 15) & 0xdb8b0000UL; /* t and c, magic vectors */ |
| y &= 0xffffffffUL; /* you may delete this line if word size = 32 */ |
| |
| /* The following line was added by Makoto Matsumoto in the 1996 |
| version to improve lower bit's correlation. Delete this line |
| to use the code published in 1994. */ |
| |
| y ^= (y >> 16); /* added to the 1994 version */ |
| |
| state->n = n + 1; |
| |
| return y; |
| } |
| |
| static double |
| tt_get_double (void * vstate) |
| { |
| return tt_get (vstate) / 4294967296.0 ; |
| } |
| |
| static void |
| tt_set (void *vstate, unsigned long int s) |
| { |
| tt_state_t *state = (tt_state_t *) vstate; |
| |
| const tt_state_t init_state = |
| {0, |
| {0x95f24dabUL, 0x0b685215UL, 0xe76ccae7UL, |
| 0xaf3ec239UL, 0x715fad23UL, 0x24a590adUL, |
| 0x69e4b5efUL, 0xbf456141UL, 0x96bc1b7bUL, |
| 0xa7bdf825UL, 0xc1de75b7UL, 0x8858a9c9UL, |
| 0x2da87693UL, 0xb657f9ddUL, 0xffdc8a9fUL, |
| 0x8121da71UL, 0x8b823ecbUL, 0x885d05f5UL, |
| 0x4e20cd47UL, 0x5a9ad5d9UL, 0x512c0c03UL, |
| 0xea857ccdUL, 0x4cc1d30fUL, 0x8891a8a1UL, |
| 0xa6b7aadbUL}}; |
| |
| |
| if (s == 0) /* default seed is given explicitly in the original code */ |
| { |
| *state = init_state; |
| } |
| else |
| { |
| int i; |
| |
| state->n = 0; |
| |
| state->x[0] = s & 0xffffffffUL; |
| |
| for (i = 1; i < N; i++) |
| state->x[i] = (69069 * state->x[i - 1]) & 0xffffffffUL; |
| } |
| |
| return; |
| } |
| |
| static const gsl_rng_type tt_type = |
| {"tt800", /* name */ |
| 0xffffffffUL, /* RAND_MAX */ |
| 0, /* RAND_MIN */ |
| sizeof (tt_state_t), |
| &tt_set, |
| &tt_get, |
| &tt_get_double}; |
| |
| const gsl_rng_type *gsl_rng_tt800 = &tt_type; |