1 /* rng/ranlxs.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 an implementation of M. Luescher's second generation
25 version of the RANLUX generator.
26
27 Thanks to Martin Luescher for providing information on this
28 generator.
29
30 */
31
32 static unsigned long int ranlxs_get (void *vstate);
33 static inline double ranlxs_get_double (void *vstate);
34 static void ranlxs_set_lux (void *state, unsigned long int s, unsigned int luxury);
35 static void ranlxs0_set (void *state, unsigned long int s);
36 static void ranlxs1_set (void *state, unsigned long int s);
37 static void ranlxs2_set (void *state, unsigned long int s);
38
39 static const int next[12] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0};
40 static const int snext[24] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
41 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0};
42
43 static const double sbase = 16777216.0; /* 2^24 */
44 static const double sone_bit = 1.0 / 16777216.0; /* 1/2^24 */
45 static const double one_bit = 1.0 / 281474976710656.0; /* 1/2^48 */
46
47 static const double shift = 268435456.0; /* 2^28 */
48
49 #define RANLUX_STEP(x1,x2,i1,i2,i3) \
50 x1=xdbl[i1] - xdbl[i2]; \
51 if (x2 < 0) \
52 { \
53 x1-=one_bit; \
54 x2+=1; \
55 } \
56 xdbl[i3]=x2
57
58 typedef struct
59 {
60 double xdbl[12], ydbl[12]; /* doubles first so they are 8-byte aligned */
61 double carry;
62 float xflt[24];
63 unsigned int ir;
64 unsigned int jr;
65 unsigned int is;
66 unsigned int is_old;
67 unsigned int pr;
68 }
69 ranlxs_state_t;
70
71 static void increment_state (ranlxs_state_t * state);
72
73 static void
increment_state(ranlxs_state_t * state)74 increment_state (ranlxs_state_t * state)
75 {
76 int k, kmax, m;
77 double x, y1, y2, y3;
78
79 float *xflt = state->xflt;
80 double *xdbl = state->xdbl;
81 double *ydbl = state->ydbl;
82 double carry = state->carry;
83 unsigned int ir = state->ir;
84 unsigned int jr = state->jr;
85
86 for (k = 0; ir > 0; ++k)
87 {
88 y1 = xdbl[jr] - xdbl[ir];
89 y2 = y1 - carry;
90 if (y2 < 0)
91 {
92 carry = one_bit;
93 y2 += 1;
94 }
95 else
96 {
97 carry = 0;
98 }
99 xdbl[ir] = y2;
100 ir = next[ir];
101 jr = next[jr];
102 }
103
104 kmax = state->pr - 12;
105
106 for (; k <= kmax; k += 12)
107 {
108 y1 = xdbl[7] - xdbl[0];
109 y1 -= carry;
110
111 RANLUX_STEP (y2, y1, 8, 1, 0);
112 RANLUX_STEP (y3, y2, 9, 2, 1);
113 RANLUX_STEP (y1, y3, 10, 3, 2);
114 RANLUX_STEP (y2, y1, 11, 4, 3);
115 RANLUX_STEP (y3, y2, 0, 5, 4);
116 RANLUX_STEP (y1, y3, 1, 6, 5);
117 RANLUX_STEP (y2, y1, 2, 7, 6);
118 RANLUX_STEP (y3, y2, 3, 8, 7);
119 RANLUX_STEP (y1, y3, 4, 9, 8);
120 RANLUX_STEP (y2, y1, 5, 10, 9);
121 RANLUX_STEP (y3, y2, 6, 11, 10);
122
123 if (y3 < 0)
124 {
125 carry = one_bit;
126 y3 += 1;
127 }
128 else
129 {
130 carry = 0;
131 }
132 xdbl[11] = y3;
133 }
134
135 kmax = state->pr;
136
137 for (; k < kmax; ++k)
138 {
139 y1 = xdbl[jr] - xdbl[ir];
140 y2 = y1 - carry;
141 if (y2 < 0)
142 {
143 carry = one_bit;
144 y2 += 1;
145 }
146 else
147 {
148 carry = 0;
149 }
150 xdbl[ir] = y2;
151 ydbl[ir] = y2 + shift;
152 ir = next[ir];
153 jr = next[jr];
154 }
155
156 ydbl[ir] = xdbl[ir] + shift;
157
158 for (k = next[ir]; k > 0;)
159 {
160 ydbl[k] = xdbl[k] + shift;
161 k = next[k];
162 }
163
164 for (k = 0, m = 0; k < 12; ++k)
165 {
166 x = xdbl[k];
167 y2 = ydbl[k] - shift;
168 if (y2 > x)
169 y2 -= sone_bit;
170 y1 = (x - y2) * sbase;
171
172 xflt[m++] = (float) y1;
173 xflt[m++] = (float) y2;
174 }
175
176 state->ir = ir;
177 state->is = 2 * ir;
178 state->is_old = 2 * ir;
179 state->jr = jr;
180 state->carry = carry;
181 }
182
183
184 static inline double
ranlxs_get_double(void * vstate)185 ranlxs_get_double (void *vstate)
186 {
187 ranlxs_state_t *state = (ranlxs_state_t *) vstate;
188
189 const unsigned int is = snext[state->is];
190
191 state->is = is;
192
193 if (is == state->is_old)
194 increment_state (state);
195
196 return state->xflt[state->is];
197 }
198
199 static unsigned long int
ranlxs_get(void * vstate)200 ranlxs_get (void *vstate)
201 {
202 return ranlxs_get_double (vstate) * 16777216.0; /* 2^24 */
203 }
204
205 static void
ranlxs_set_lux(void * vstate,unsigned long int s,unsigned int luxury)206 ranlxs_set_lux (void *vstate, unsigned long int s, unsigned int luxury)
207 {
208 ranlxs_state_t *state = (ranlxs_state_t *) vstate;
209
210 int ibit, jbit, i, k, m, xbit[31];
211 double x, y;
212
213 long int seed;
214
215 if (s == 0)
216 s = 1; /* default seed is 1 */
217
218 seed = s;
219
220 i = seed & 0xFFFFFFFFUL;
221
222 for (k = 0; k < 31; ++k)
223 {
224 xbit[k] = i % 2;
225 i /= 2;
226 }
227
228 ibit = 0;
229 jbit = 18;
230
231 for (k = 0; k < 12; ++k)
232 {
233 x = 0;
234
235 for (m = 1; m <= 48; ++m)
236 {
237 y = (double) xbit[ibit];
238 x += x + y;
239 xbit[ibit] = (xbit[ibit] + xbit[jbit]) % 2;
240 ibit = (ibit + 1) % 31;
241 jbit = (jbit + 1) % 31;
242 }
243 state->xdbl[k] = one_bit * x;
244 }
245
246 state->carry = 0;
247 state->ir = 0;
248 state->jr = 7;
249 state->is = 23;
250 state->is_old = 0;
251 state->pr = luxury;
252 }
253
254 static void
ranlxs0_set(void * vstate,unsigned long int s)255 ranlxs0_set (void *vstate, unsigned long int s)
256 {
257 ranlxs_set_lux (vstate, s, 109);
258 }
259
260 void
ranlxs1_set(void * vstate,unsigned long int s)261 ranlxs1_set (void *vstate, unsigned long int s)
262 {
263 ranlxs_set_lux (vstate, s, 202);
264 }
265
266 static void
ranlxs2_set(void * vstate,unsigned long int s)267 ranlxs2_set (void *vstate, unsigned long int s)
268 {
269 ranlxs_set_lux (vstate, s, 397);
270 }
271
272
273 static const gsl_rng_type ranlxs0_type =
274 {"ranlxs0", /* name */
275 0x00ffffffUL, /* RAND_MAX */
276 0, /* RAND_MIN */
277 sizeof (ranlxs_state_t),
278 &ranlxs0_set,
279 &ranlxs_get,
280 &ranlxs_get_double};
281
282 static const gsl_rng_type ranlxs1_type =
283 {"ranlxs1", /* name */
284 0x00ffffffUL, /* RAND_MAX */
285 0, /* RAND_MIN */
286 sizeof (ranlxs_state_t),
287 &ranlxs1_set,
288 &ranlxs_get,
289 &ranlxs_get_double};
290
291 static const gsl_rng_type ranlxs2_type =
292 {"ranlxs2", /* name */
293 0x00ffffffUL, /* RAND_MAX */
294 0, /* RAND_MIN */
295 sizeof (ranlxs_state_t),
296 &ranlxs2_set,
297 &ranlxs_get,
298 &ranlxs_get_double};
299
300 const gsl_rng_type *gsl_rng_ranlxs0 = &ranlxs0_type;
301 const gsl_rng_type *gsl_rng_ranlxs1 = &ranlxs1_type;
302 const gsl_rng_type *gsl_rng_ranlxs2 = &ranlxs2_type;
303