3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
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.
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.
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.
22 #include <gsl/gsl_rng.h>
24 /* This is an implementation of M. Luescher's second generation
25 version of the RANLUX generator.
27 Thanks to Martin Luescher for providing information on this
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);
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};
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 */
47 static const double shift = 268435456.0; /* 2^28 */
49 #define RANLUX_STEP(x1,x2,i1,i2,i3) \
50 x1=xdbl[i1] - xdbl[i2]; \
60 double xdbl[12], ydbl[12]; /* doubles first so they are 8-byte aligned */
71 static void increment_state (ranlxs_state_t * state);
74 increment_state (ranlxs_state_t * state)
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;
86 for (k = 0; ir > 0; ++k)
88 y1 = xdbl[jr] - xdbl[ir];
104 kmax = state->pr - 12;
106 for (; k <= kmax; k += 12)
108 y1 = xdbl[7] - xdbl[0];
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);
137 for (; k < kmax; ++k)
139 y1 = xdbl[jr] - xdbl[ir];
151 ydbl[ir] = y2 + shift;
156 ydbl[ir] = xdbl[ir] + shift;
158 for (k = next[ir]; k > 0;)
160 ydbl[k] = xdbl[k] + shift;
164 for (k = 0, m = 0; k < 12; ++k)
167 y2 = ydbl[k] - shift;
170 y1 = (x - y2) * sbase;
172 xflt[m++] = (float) y1;
173 xflt[m++] = (float) y2;
178 state->is_old = 2 * ir;
180 state->carry = carry;
185 ranlxs_get_double (void *vstate)
187 ranlxs_state_t *state = (ranlxs_state_t *) vstate;
189 const unsigned int is = snext[state->is];
193 if (is == state->is_old)
194 increment_state (state);
196 return state->xflt[state->is];
199 static unsigned long int
200 ranlxs_get (void *vstate)
202 return ranlxs_get_double (vstate) * 16777216.0; /* 2^24 */
206 ranlxs_set_lux (void *vstate, unsigned long int s, unsigned int luxury)
208 ranlxs_state_t *state = (ranlxs_state_t *) vstate;
210 int ibit, jbit, i, k, m, xbit[31];
216 s = 1; /* default seed is 1 */
220 i = seed & 0xFFFFFFFFUL;
222 for (k = 0; k < 31; ++k)
231 for (k = 0; k < 12; ++k)
235 for (m = 1; m <= 48; ++m)
237 y = (double) xbit[ibit];
239 xbit[ibit] = (xbit[ibit] + xbit[jbit]) % 2;
240 ibit = (ibit + 1) % 31;
241 jbit = (jbit + 1) % 31;
243 state->xdbl[k] = one_bit * x;
255 ranlxs0_set (void *vstate, unsigned long int s)
257 ranlxs_set_lux (vstate, s, 109);
261 ranlxs1_set (void *vstate, unsigned long int s)
263 ranlxs_set_lux (vstate, s, 202);
267 ranlxs2_set (void *vstate, unsigned long int s)
269 ranlxs_set_lux (vstate, s, 397);
273 static const gsl_rng_type ranlxs0_type =
274 {"ranlxs0", /* name */
275 0x00ffffffUL, /* RAND_MAX */
277 sizeof (ranlxs_state_t),
282 static const gsl_rng_type ranlxs1_type =
283 {"ranlxs1", /* name */
284 0x00ffffffUL, /* RAND_MAX */
286 sizeof (ranlxs_state_t),
291 static const gsl_rng_type ranlxs2_type =
292 {"ranlxs2", /* name */
293 0x00ffffffUL, /* RAND_MAX */
295 sizeof (ranlxs_state_t),
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;