|
Packit |
67cb25 |
/* This program is free software; you can redistribute it and/or
|
|
Packit |
67cb25 |
modify it under the terms of the GNU General Public License as
|
|
Packit |
67cb25 |
published by the Free Software Foundation; either version 3 of the
|
|
Packit |
67cb25 |
License, or (at your option) any later version.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This program is distributed in the hope that it will be useful,
|
|
Packit |
67cb25 |
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
General Public License for more details. You should have received
|
|
Packit |
67cb25 |
a copy of the GNU General Public License along with this library; if
|
|
Packit |
67cb25 |
not, write to the Free Software Foundation, Inc., 51 Franklin Street,
|
|
Packit |
67cb25 |
Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
From Robert M. Ziff, "Four-tap shift-register-sequence
|
|
Packit |
67cb25 |
random-number generators," Computers in Physics 12(4), Jul/Aug
|
|
Packit |
67cb25 |
1998, pp 385-392. A generalized feedback shift-register (GFSR)
|
|
Packit |
67cb25 |
is basically an xor-sum of particular past lagged values. A
|
|
Packit |
67cb25 |
four-tap register looks like:
|
|
Packit |
67cb25 |
ra[nd] = ra[nd-A] ^ ra[nd-B] ^ ra[nd-C] ^ ra[nd-D]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Ziff notes that "it is now widely known" that two-tap registers
|
|
Packit |
67cb25 |
have serious flaws, the most obvious one being the three-point
|
|
Packit |
67cb25 |
correlation that comes from the defn of the generator. Nice
|
|
Packit |
67cb25 |
mathematical properties can be derived for GFSR's, and numerics
|
|
Packit |
67cb25 |
bears out the claim that 4-tap GFSR's with appropriately chosen
|
|
Packit |
67cb25 |
offsets are as random as can be measured, using the author's test.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This implementation uses the values suggested the the author's
|
|
Packit |
67cb25 |
example on p392, but altered to fit the GSL framework. The "state"
|
|
Packit |
67cb25 |
is 2^14 longs, or 64Kbytes; 2^14 is the smallest power of two that
|
|
Packit |
67cb25 |
is larger than D, the largest offset. We really only need a state
|
|
Packit |
67cb25 |
with the last D values, but by going to a power of two, we can do a
|
|
Packit |
67cb25 |
masking operation instead of a modulo, and this is presumably
|
|
Packit |
67cb25 |
faster, though I haven't actually tried it. The article actually
|
|
Packit |
67cb25 |
suggested a short/fast hack:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define RandomInteger (++nd, ra[nd&M]=ra[(nd-A)&M]\
|
|
Packit |
67cb25 |
^ra[(nd-B)&M]^ra[(nd-C)&M]^ra[(nd-D)&M])
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
so that (as long as you've defined nd,ra[M+1]), then you ca do things
|
|
Packit |
67cb25 |
like: 'if (RandomInteger < p) {...}'.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Note that n&M varies from 0 to M, *including* M, so that the
|
|
Packit |
67cb25 |
array has to be of size M+1. Since M+1 is a power of two, n&M
|
|
Packit |
67cb25 |
is a potentially quicker implementation of the equivalent n%(M+1).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This implementation copyright (C) 1998 James Theiler, based on
|
|
Packit |
67cb25 |
the example mt.c in the GSL, as implemented by Brian Gough.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_rng.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline unsigned long int gfsr4_get (void *vstate);
|
|
Packit |
67cb25 |
static double gfsr4_get_double (void *vstate);
|
|
Packit |
67cb25 |
static void gfsr4_set (void *state, unsigned long int s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Magic numbers */
|
|
Packit |
67cb25 |
#define A 471
|
|
Packit |
67cb25 |
#define B 1586
|
|
Packit |
67cb25 |
#define C 6988
|
|
Packit |
67cb25 |
#define D 9689
|
|
Packit |
67cb25 |
#define M 16383 /* = 2^14-1 */
|
|
Packit |
67cb25 |
/* #define M 0x0003fff */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
typedef struct
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int nd;
|
|
Packit |
67cb25 |
unsigned long ra[M+1];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
gfsr4_state_t;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline unsigned long
|
|
Packit |
67cb25 |
gfsr4_get (void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gfsr4_state_t *state = (gfsr4_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->nd = ((state->nd)+1)&M;
|
|
Packit |
67cb25 |
return state->ra[(state->nd)] =
|
|
Packit |
67cb25 |
state->ra[((state->nd)+(M+1-A))&M]^
|
|
Packit |
67cb25 |
state->ra[((state->nd)+(M+1-B))&M]^
|
|
Packit |
67cb25 |
state->ra[((state->nd)+(M+1-C))&M]^
|
|
Packit |
67cb25 |
state->ra[((state->nd)+(M+1-D))&M];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
gfsr4_get_double (void * vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return gfsr4_get (vstate) / 4294967296.0 ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
gfsr4_set (void *vstate, unsigned long int s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gfsr4_state_t *state = (gfsr4_state_t *) vstate;
|
|
Packit |
67cb25 |
int i, j;
|
|
Packit |
67cb25 |
/* Masks for turning on the diagonal bit and turning off the
|
|
Packit |
67cb25 |
leftmost bits */
|
|
Packit |
67cb25 |
unsigned long int msb = 0x80000000UL;
|
|
Packit |
67cb25 |
unsigned long int mask = 0xffffffffUL;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s == 0)
|
|
Packit |
67cb25 |
s = 4357; /* the default seed is 4357 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* We use the congruence s_{n+1} = (69069*s_n) mod 2^32 to
|
|
Packit |
67cb25 |
initialize the state. This works because ANSI-C unsigned long
|
|
Packit |
67cb25 |
integer arithmetic is automatically modulo 2^32 (or a higher
|
|
Packit |
67cb25 |
power of two), so we can safely ignore overflow. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define LCG(n) ((69069 * n) & 0xffffffffUL)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Brian Gough suggests this to avoid low-order bit correlations */
|
|
Packit |
67cb25 |
for (i = 0; i <= M; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
unsigned long t = 0 ;
|
|
Packit |
67cb25 |
unsigned long bit = msb ;
|
|
Packit |
67cb25 |
for (j = 0; j < 32; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
s = LCG(s) ;
|
|
Packit |
67cb25 |
if (s & msb)
|
|
Packit |
67cb25 |
t |= bit ;
|
|
Packit |
67cb25 |
bit >>= 1 ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
state->ra[i] = t ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Perform the "orthogonalization" of the matrix */
|
|
Packit |
67cb25 |
/* Based on the orthogonalization used in r250, as suggested initially
|
|
Packit |
67cb25 |
* by Kirkpatrick and Stoll, and pointed out to me by Brian Gough
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* BJG: note that this orthogonalisation doesn't have any effect
|
|
Packit |
67cb25 |
here because the the initial 6695 elements do not participate in
|
|
Packit |
67cb25 |
the calculation. For practical purposes this orthogonalisation
|
|
Packit |
67cb25 |
is somewhat irrelevant, because the probability of the original
|
|
Packit |
67cb25 |
sequence being degenerate should be exponentially small. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i=0; i<32; ++i) {
|
|
Packit |
67cb25 |
int k=7+i*3;
|
|
Packit |
67cb25 |
state->ra[k] &= mask; /* Turn off bits left of the diagonal */
|
|
Packit |
67cb25 |
state->ra[k] |= msb; /* Turn on the diagonal bit */
|
|
Packit |
67cb25 |
mask >>= 1;
|
|
Packit |
67cb25 |
msb >>= 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->nd = i;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static const gsl_rng_type gfsr4_type =
|
|
Packit |
67cb25 |
{"gfsr4", /* name */
|
|
Packit |
67cb25 |
0xffffffffUL, /* RAND_MAX */
|
|
Packit |
67cb25 |
0, /* RAND_MIN */
|
|
Packit |
67cb25 |
sizeof (gfsr4_state_t),
|
|
Packit |
67cb25 |
&gfsr4_set,
|
|
Packit |
67cb25 |
&gfsr4_get,
|
|
Packit |
67cb25 |
&gfsr4_get_double};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const gsl_rng_type *gsl_rng_gfsr4 = &gfsr4_type;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|