| /* 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; |