diff options
Diffstat (limited to 'gsl-1.9/rng/knuthran2002.c')
| -rw-r--r-- | gsl-1.9/rng/knuthran2002.c | 189 |
1 files changed, 189 insertions, 0 deletions
diff --git a/gsl-1.9/rng/knuthran2002.c b/gsl-1.9/rng/knuthran2002.c new file mode 100644 index 0000000..838b351 --- /dev/null +++ b/gsl-1.9/rng/knuthran2002.c @@ -0,0 +1,189 @@ +/* rng/knuthran2002.c + * + * Copyright (C) 2007 Brian Gough + * 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, Section 3.6 + * Third Edition, Addison-Wesley, + * + * The modifications introduced in the 9th printing (2002) are + * included here; there's no backwards compatibility with the + * original. [ see http://www-cs-faculty.stanford.edu/~knuth/taocp.html ] + * + */ + +#include <config.h> +#include <stdlib.h> +#include <gsl/gsl_rng.h> + +#define BUFLEN 1009 /* 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 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 (long int aa[], unsigned int n, + 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; + long int aa[BUFLEN]; + long int ran_x[KK]; /* the generator state */ +} +ran_state_t; + +static inline void +ran_array (long int aa[], unsigned int n, 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; + unsigned long int v; + + if (i == 0) + { + /* fill buffer with new random numbers */ + ran_array (state->aa, BUFLEN, state->ran_x); + } + + v = state->aa[i]; + + state->i = (i + 1) % KK; + + return v; +} + +static double +ran_get_double (void *vstate) +{ + ran_state_t *state = (ran_state_t *) vstate; + + return ran_get (state) / 1073741824.0; /* 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; + + if (s == 0 ) + s = 314159; /* default seed used by Knuth */ + + ss = (s + 2)&(MM-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; + } + 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]; + x[j + j - 1] = 0; + } + + for (j = KK + KK - 2; j >= KK; 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 */ + x[LL] = mod_diff (x[LL], x[KK]); + } + + if (ss) + ss >>= 1; + else + t--; + } + + 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]; + + + for (j = 0; j< 10; j++) + ran_array(x, KK+KK-1, state->ran_x); /* warm things up */ + + state->i = 0; + + return; +} + +static const gsl_rng_type ran_type = { + "knuthran2002", /* name */ + 0x3fffffffUL, /* RAND_MAX = (2 ^ 30) - 1 */ + 0, /* RAND_MIN */ + sizeof (ran_state_t), + &ran_set, + &ran_get, + &ran_get_double +}; + +const gsl_rng_type *gsl_rng_knuthran2002 = &ran_type; |