1 /* rng/tt.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include "gsl__config.h"
21 #include <stdlib.h>
22 #include "gsl_rng.h"
23
24 /* This is the TT800 twisted GSFR generator for 32 bit integers. It
25 has been superceded by MT19937 (mt.c). The period is 2^800.
26
27 This implementation is based on tt800.c, July 8th 1996 version by
28 M. Matsumoto, email: matumoto@math.keio.ac.jp
29
30 From: Makoto Matsumoto and Yoshiharu Kurita, "Twisted GFSR
31 Generators II", ACM Transactions on Modelling and Computer
32 Simulation, Vol. 4, No. 3, 1994, pages 254-266. */
33
34 static inline unsigned long int tt_get (void *vstate);
35 static double tt_get_double (void *vstate);
36 static void tt_set (void *state, unsigned long int s);
37
38 #define N 25
39 #define M 7
40
41 typedef struct
42 {
43 int n;
44 unsigned long int x[N];
45 }
46 tt_state_t;
47
48 static inline unsigned long int
tt_get(void * vstate)49 tt_get (void *vstate)
50 {
51 tt_state_t *state = (tt_state_t *) vstate;
52
53 /* this is the magic vector, a */
54
55 const unsigned long mag01[2] =
56 {0x00000000, 0x8ebfd028UL};
57 unsigned long int y;
58 unsigned long int *const x = state->x;
59 int n = state->n;
60
61 if (n >= N)
62 {
63 int i;
64 for (i = 0; i < N - M; i++)
65 {
66 x[i] = x[i + M] ^ (x[i] >> 1) ^ mag01[x[i] % 2];
67 }
68 for (; i < N; i++)
69 {
70 x[i] = x[i + (M - N)] ^ (x[i] >> 1) ^ mag01[x[i] % 2];
71 };
72 n = 0;
73 }
74
75 y = x[n];
76 y ^= (y << 7) & 0x2b5b2500UL; /* s and b, magic vectors */
77 y ^= (y << 15) & 0xdb8b0000UL; /* t and c, magic vectors */
78 y &= 0xffffffffUL; /* you may delete this line if word size = 32 */
79
80 /* The following line was added by Makoto Matsumoto in the 1996
81 version to improve lower bit's correlation. Delete this line
82 to use the code published in 1994. */
83
84 y ^= (y >> 16); /* added to the 1994 version */
85
86 state->n = n + 1;
87
88 return y;
89 }
90
91 static double
tt_get_double(void * vstate)92 tt_get_double (void * vstate)
93 {
94 return tt_get (vstate) / 4294967296.0 ;
95 }
96
97 static void
tt_set(void * vstate,unsigned long int s)98 tt_set (void *vstate, unsigned long int s)
99 {
100 tt_state_t *state = (tt_state_t *) vstate;
101
102 const tt_state_t init_state =
103 {0,
104 {0x95f24dabUL, 0x0b685215UL, 0xe76ccae7UL,
105 0xaf3ec239UL, 0x715fad23UL, 0x24a590adUL,
106 0x69e4b5efUL, 0xbf456141UL, 0x96bc1b7bUL,
107 0xa7bdf825UL, 0xc1de75b7UL, 0x8858a9c9UL,
108 0x2da87693UL, 0xb657f9ddUL, 0xffdc8a9fUL,
109 0x8121da71UL, 0x8b823ecbUL, 0x885d05f5UL,
110 0x4e20cd47UL, 0x5a9ad5d9UL, 0x512c0c03UL,
111 0xea857ccdUL, 0x4cc1d30fUL, 0x8891a8a1UL,
112 0xa6b7aadbUL}};
113
114
115 if (s == 0) /* default seed is given explicitly in the original code */
116 {
117 *state = init_state;
118 }
119 else
120 {
121 int i;
122
123 state->n = 0;
124
125 state->x[0] = s & 0xffffffffUL;
126
127 for (i = 1; i < N; i++)
128 state->x[i] = (69069 * state->x[i - 1]) & 0xffffffffUL;
129 }
130
131 return;
132 }
133
134 static const gsl_rng_type tt_type =
135 {"tt800", /* name */
136 0xffffffffUL, /* RAND_MAX */
137 0, /* RAND_MIN */
138 sizeof (tt_state_t),
139 &tt_set,
140 &tt_get,
141 &tt_get_double};
142
143 const gsl_rng_type *gsl_rng_tt800 = &tt_type;
144