diff options
Diffstat (limited to 'gsl-1.9/rng/knuthran.c')
| -rw-r--r-- | gsl-1.9/rng/knuthran.c | 177 |
1 files changed, 177 insertions, 0 deletions
diff --git a/gsl-1.9/rng/knuthran.c b/gsl-1.9/rng/knuthran.c new file mode 100644 index 0000000..e83678b --- /dev/null +++ b/gsl-1.9/rng/knuthran.c @@ -0,0 +1,177 @@ +/* rng/knuthran.c + * + * Copyright (C) 2001 Brian Gough, Carlo Perassi + * + * 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. + */ + +/* + * This generator is taken from + * + * Donald E. Knuth + * The Art of Computer Programming + * Volume 2 + * Third Edition + * Addison-Wesley + * Section 3.6 + * + * The comments are taken from the book + * Our comments are signed + */ + +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_rng.h> + +#define BUFLEN 2009 /* [Brian]: length of the buffer aa[] */ +#define KK 100 /* the long lag */ +#define LL 37 /* the short lag */ +#define MM (1L << 30) /* the modulus */ +#define TT 70 /* guaranteed separation between streams */ + +#define evenize(x) ((x) & (MM - 2)) /* make x even */ +#define is_odd(x) ((x) & 1) /* the units bit of x */ +#define mod_diff(x, y) (((x) - (y)) & (MM - 1)) /* (x - y) mod MM */ + +static inline void ran_array (unsigned long int aa[], unsigned int n, + unsigned long int ran_x[]); +static inline unsigned long int ran_get (void *vstate); +static double ran_get_double (void *vstate); +static void ran_set (void *state, unsigned long int s); + +typedef struct +{ + unsigned int i; + unsigned long int aa[BUFLEN]; /* [Carlo]: I can't pass n to ran_array like + Knuth does */ + unsigned long int ran_x[KK]; /* the generator state */ +} +ran_state_t; + +static inline void +ran_array (unsigned long int aa[], unsigned int n, unsigned long int ran_x[]) +{ + unsigned int i; + unsigned int j; + + for (j = 0; j < KK; j++) + aa[j] = ran_x[j]; + + for (; j < n; j++) + aa[j] = mod_diff (aa[j - KK], aa[j - LL]); + + for (i = 0; i < LL; i++, j++) + ran_x[i] = mod_diff (aa[j - KK], aa[j - LL]); + + for (; i < KK; i++, j++) + ran_x[i] = mod_diff (aa[j - KK], ran_x[i - LL]); +} + +static inline unsigned long int +ran_get (void *vstate) +{ + ran_state_t *state = (ran_state_t *) vstate; + + unsigned int i = state->i; + + if (i == 0) + { + /* fill buffer with new random numbers */ + ran_array (state->aa, BUFLEN, state->ran_x); + } + + state->i = (i + 1) % BUFLEN; + + return state->aa[i]; +} + +static double +ran_get_double (void *vstate) +{ + ran_state_t *state = (ran_state_t *) vstate; + + return ran_get (state) / 1073741824.0; /* [Carlo]: RAND_MAX + 1 */ +} + +static void +ran_set (void *vstate, unsigned long int s) +{ + ran_state_t *state = (ran_state_t *) vstate; + + long x[KK + KK - 1]; /* the preparation buffer */ + + register int j; + register int t; + register long ss = evenize (s + 2); + + for (j = 0; j < KK; j++) + { + x[j] = ss; /* bootstrap the buffer */ + ss <<= 1; + if (ss >= MM) /* cyclic shift 29 bits */ + ss -= MM - 2; + } + for (; j < KK + KK - 1; j++) + x[j] = 0; + x[1]++; /* make x[1] (and only x[1]) odd */ + ss = s & (MM - 1); + t = TT - 1; + while (t) + { + for (j = KK - 1; j > 0; j--) /* square */ + x[j + j] = x[j]; + for (j = KK + KK - 2; j > KK - LL; j -= 2) + x[KK + KK - 1 - j] = evenize (x[j]); + for (j = KK + KK - 2; j >= KK; j--) + if (is_odd (x[j])) + { + x[j - (KK - LL)] = mod_diff (x[j - (KK - LL)], x[j]); + x[j - KK] = mod_diff (x[j - KK], x[j]); + } + if (is_odd (ss)) + { /* multiply by "z" */ + for (j = KK; j > 0; j--) + x[j] = x[j - 1]; + x[0] = x[KK]; /* shift the buffer cyclically */ + if (is_odd (x[KK])) + x[LL] = mod_diff (x[LL], x[KK]); + } + if (ss) + ss >>= 1; + else + t--; + } + + state->i = 0; + + for (j = 0; j < LL; j++) + state->ran_x[j + KK - LL] = x[j]; + for (; j < KK; j++) + state->ran_x[j - LL] = x[j]; + + return; +} + +static const gsl_rng_type ran_type = { + "knuthran", /* name */ + 0x3fffffffUL, /* RAND_MAX *//* [Carlo]: (2 ^ 30) - 1 */ + 0, /* RAND_MIN */ + sizeof (ran_state_t), + &ran_set, + &ran_get, + &ran_get_double +}; + +const gsl_rng_type *gsl_rng_knuthran = &ran_type; |