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