diff options
Diffstat (limited to 'gsl-1.9/rng/r250.c')
| -rw-r--r-- | gsl-1.9/rng/r250.c | 172 |
1 files changed, 172 insertions, 0 deletions
diff --git a/gsl-1.9/rng/r250.c b/gsl-1.9/rng/r250.c new file mode 100644 index 0000000..e5b4df1 --- /dev/null +++ b/gsl-1.9/rng/r250.c @@ -0,0 +1,172 @@ +/* rng/r250.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 a shift-register random number generator. The sequence is + + x_n = x_{n-103} ^ x_{n-250} ("^" means XOR) + + defined on 32-bit words. + + BJG: Note that this implementation actually uses the sequence, x_n + = x_{n-147} ^ x_{n-250} which generates the outputs in + time-reversed order but is otherwise completely equivalent. + + The first 250 elements x_1 .. x_250 are first initialized as x_n = + s_n, where s_n = (69069*s_{n-1}) mod 2^32 and s_0=s is the + user-supplied seed. To ensure that the sequence does not lie on a + subspace we force 32 of the entries to be linearly independent. We + take the 32 elements x[3], x[10], x[17], x[24], ..., 213 and apply + the following operations, + + x[3] &= 11111111111111111111111111111111 + x[3] |= 10000000000000000000000000000000 + x[10] &= 01111111111111111111111111111111 + x[10] |= 01000000000000000000000000000000 + x[17] &= 00111111111111111111111111111111 + x[17] |= 00100000000000000000000000000000 + .... ... + x[206] &= 00000000000000000000000000000111 + x[206] |= 00000000000000000000000000000100 + x[213] &= 00000000000000000000000000000011 + x[213] |= 00000000000000000000000000000010 + x[220] &= 00000000000000000000000000000001 + x[220] |= 00000000000000000000000000000001 + + i.e. if we consider the bits of the 32 elements as forming a 32x32 + array then we are setting the diagonal bits of the array to one and + masking the lower triangle below the diagonal to zero. + + With this initialization procedure the theoretical value of + x_{10001} is 1100653588 for s = 1 (Actually I got this by running + the original code). The subscript 10001 means (1) seed the + generator with s = 1 and then do 10000 actual iterations. + + The period of this generator is about 2^250. + + The algorithm works for any number of bits. It is implemented here + for 32 bits. + + From: S. Kirkpatrick and E. Stoll, "A very fast shift-register + sequence random number generator", Journal of Computational Physics, + 40, 517-526 (1981). */ + +static inline unsigned long int r250_get (void *vstate); +static double r250_get_double (void *vstate); +static void r250_set (void *state, unsigned long int s); + +typedef struct + { + int i; + unsigned long x[250]; + } +r250_state_t; + +static inline unsigned long int +r250_get (void *vstate) +{ + r250_state_t *state = (r250_state_t *) vstate; + unsigned long int k; + int j; + + int i = state->i; + + if (i >= 147) + { + j = i - 147; + } + else + { + j = i + 103; + } + + k = state->x[i] ^ state->x[j]; + state->x[i] = k; + + if (i >= 249) + { + state->i = 0; + } + else + { + state->i = i + 1; + } + + return k; +} + +static double +r250_get_double (void *vstate) +{ + return r250_get (vstate) / 4294967296.0 ; +} + +static void +r250_set (void *vstate, unsigned long int s) +{ + r250_state_t *state = (r250_state_t *) vstate; + + int i; + + if (s == 0) + s = 1; /* default seed is 1 */ + + state->i = 0; + +#define LCG(n) ((69069 * n) & 0xffffffffUL) + + for (i = 0; i < 250; i++) /* Fill the buffer */ + { + s = LCG (s); + state->x[i] = s; + } + + { + /* Masks for turning on the diagonal bit and turning off the + leftmost bits */ + + unsigned long int msb = 0x80000000UL; + unsigned long int mask = 0xffffffffUL; + + for (i = 0; i < 32; i++) + { + int k = 7 * i + 3; /* Select a word to operate on */ + state->x[k] &= mask; /* Turn off bits left of the diagonal */ + state->x[k] |= msb; /* Turn on the diagonal bit */ + mask >>= 1; + msb >>= 1; + } + } + + return; +} + +static const gsl_rng_type r250_type = +{"r250", /* name */ + 0xffffffffUL, /* RAND_MAX */ + 0, /* RAND_MIN */ + sizeof (r250_state_t), + &r250_set, + &r250_get, + &r250_get_double}; + +const gsl_rng_type *gsl_rng_r250 = &r250_type; |