1 /* rng/ranlxd.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 <config.h>
21 #include <stdlib.h>
22 #include <gsl/gsl_rng.h>
23
24 /* This is an implementation of Martin Luescher's second generation
25 double-precision (48-bit) version of the RANLUX generator.
26
27 Thanks to Martin Luescher for providing information on this
28 generator.
29
30 */
31
32 static inline unsigned long int ranlxd_get (void *vstate);
33 static double ranlxd_get_double (void *vstate);
34 static void ranlxd_set_lux (void *state, unsigned long int s, unsigned int luxury);
35 static void ranlxd1_set (void *state, unsigned long int s);
36 static void ranlxd2_set (void *state, unsigned long int s);
37
38 static const int next[12] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0};
39
40 static const double one_bit = 1.0 / 281474976710656.0; /* 1/2^48 */
41
42 #define RANLUX_STEP(x1,x2,i1,i2,i3) \
43 x1=xdbl[i1] - xdbl[i2]; \
44 if (x2 < 0) \
45 { \
46 x1-=one_bit; \
47 x2+=1; \
48 } \
49 xdbl[i3]=x2
50
51 typedef struct
52 {
53 double xdbl[12];
54 double carry;
55 unsigned int ir;
56 unsigned int jr;
57 unsigned int ir_old;
58 unsigned int pr;
59 }
60 ranlxd_state_t;
61
62 static inline void increment_state (ranlxd_state_t * state);
63
64 static inline void
increment_state(ranlxd_state_t * state)65 increment_state (ranlxd_state_t * state)
66 {
67 int k, kmax;
68 double y1, y2, y3;
69
70 double *xdbl = state->xdbl;
71 double carry = state->carry;
72 unsigned int ir = state->ir;
73 unsigned int jr = state->jr;
74
75 for (k = 0; ir > 0; ++k)
76 {
77 y1 = xdbl[jr] - xdbl[ir];
78 y2 = y1 - carry;
79 if (y2 < 0)
80 {
81 carry = one_bit;
82 y2 += 1;
83 }
84 else
85 {
86 carry = 0;
87 }
88 xdbl[ir] = y2;
89 ir = next[ir];
90 jr = next[jr];
91 }
92
93 kmax = state->pr - 12;
94
95 for (; k <= kmax; k += 12)
96 {
97 y1 = xdbl[7] - xdbl[0];
98 y1 -= carry;
99
100 RANLUX_STEP (y2, y1, 8, 1, 0);
101 RANLUX_STEP (y3, y2, 9, 2, 1);
102 RANLUX_STEP (y1, y3, 10, 3, 2);
103 RANLUX_STEP (y2, y1, 11, 4, 3);
104 RANLUX_STEP (y3, y2, 0, 5, 4);
105 RANLUX_STEP (y1, y3, 1, 6, 5);
106 RANLUX_STEP (y2, y1, 2, 7, 6);
107 RANLUX_STEP (y3, y2, 3, 8, 7);
108 RANLUX_STEP (y1, y3, 4, 9, 8);
109 RANLUX_STEP (y2, y1, 5, 10, 9);
110 RANLUX_STEP (y3, y2, 6, 11, 10);
111
112 if (y3 < 0)
113 {
114 carry = one_bit;
115 y3 += 1;
116 }
117 else
118 {
119 carry = 0;
120 }
121 xdbl[11] = y3;
122 }
123
124 kmax = state->pr;
125
126 for (; k < kmax; ++k)
127 {
128 y1 = xdbl[jr] - xdbl[ir];
129 y2 = y1 - carry;
130 if (y2 < 0)
131 {
132 carry = one_bit;
133 y2 += 1;
134 }
135 else
136 {
137 carry = 0;
138 }
139 xdbl[ir] = y2;
140 ir = next[ir];
141 jr = next[jr];
142 }
143 state->ir = ir;
144 state->ir_old = ir;
145 state->jr = jr;
146 state->carry = carry;
147 }
148
149 static inline unsigned long int
ranlxd_get(void * vstate)150 ranlxd_get (void *vstate)
151 {
152 return ranlxd_get_double (vstate) * 4294967296.0; /* 2^32 */
153 }
154
155 static double
ranlxd_get_double(void * vstate)156 ranlxd_get_double (void *vstate)
157 {
158 ranlxd_state_t *state = (ranlxd_state_t *) vstate;
159
160 int ir = state->ir;
161
162 state->ir = next[ir];
163
164 if (state->ir == state->ir_old)
165 increment_state (state);
166
167 return state->xdbl[state->ir];
168 }
169
170 static void
ranlxd_set_lux(void * vstate,unsigned long int s,unsigned int luxury)171 ranlxd_set_lux (void *vstate, unsigned long int s, unsigned int luxury)
172 {
173 ranlxd_state_t *state = (ranlxd_state_t *) vstate;
174
175 int ibit, jbit, i, k, l, xbit[31];
176 double x, y;
177
178 long int seed;
179
180 if (s == 0)
181 s = 1; /* default seed is 1 */
182
183 seed = s;
184
185 i = seed & 0xFFFFFFFFUL;
186
187 for (k = 0; k < 31; ++k)
188 {
189 xbit[k] = i % 2;
190 i /= 2;
191 }
192
193 ibit = 0;
194 jbit = 18;
195
196 for (k = 0; k < 12; ++k)
197 {
198 x = 0;
199
200 for (l = 1; l <= 48; ++l)
201 {
202 y = (double) ((xbit[ibit] + 1) % 2);
203 x += x + y;
204 xbit[ibit] = (xbit[ibit] + xbit[jbit]) % 2;
205 ibit = (ibit + 1) % 31;
206 jbit = (jbit + 1) % 31;
207 }
208 state->xdbl[k] = one_bit * x;
209 }
210
211 state->carry = 0;
212 state->ir = 11;
213 state->jr = 7;
214 state->ir_old = 0;
215 state->pr = luxury;
216 }
217
218 static void
ranlxd1_set(void * vstate,unsigned long int s)219 ranlxd1_set (void *vstate, unsigned long int s)
220 {
221 ranlxd_set_lux (vstate, s, 202);
222 }
223
224 static void
ranlxd2_set(void * vstate,unsigned long int s)225 ranlxd2_set (void *vstate, unsigned long int s)
226 {
227 ranlxd_set_lux (vstate, s, 397);
228 }
229
230 static const gsl_rng_type ranlxd1_type =
231 {"ranlxd1", /* name */
232 0xffffffffUL, /* RAND_MAX */
233 0, /* RAND_MIN */
234 sizeof (ranlxd_state_t),
235 &ranlxd1_set,
236 &ranlxd_get,
237 &ranlxd_get_double};
238
239 static const gsl_rng_type ranlxd2_type =
240 {"ranlxd2", /* name */
241 0xffffffffUL, /* RAND_MAX */
242 0, /* RAND_MIN */
243 sizeof (ranlxd_state_t),
244 &ranlxd2_set,
245 &ranlxd_get,
246 &ranlxd_get_double};
247
248 const gsl_rng_type *gsl_rng_ranlxd1 = &ranlxd1_type;
249 const gsl_rng_type *gsl_rng_ranlxd2 = &ranlxd2_type;
250