|
Packit |
67cb25 |
/* rng/r250.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* 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.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
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 |
/* This is a shift-register random number generator. The sequence is
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x_n = x_{n-103} ^ x_{n-250} ("^" means XOR)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
defined on 32-bit words.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
BJG: Note that this implementation actually uses the sequence, x_n
|
|
Packit |
67cb25 |
= x_{n-147} ^ x_{n-250} which generates the outputs in
|
|
Packit |
67cb25 |
time-reversed order but is otherwise completely equivalent.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The first 250 elements x_1 .. x_250 are first initialized as x_n =
|
|
Packit |
67cb25 |
s_n, where s_n = (69069*s_{n-1}) mod 2^32 and s_0=s is the
|
|
Packit |
67cb25 |
user-supplied seed. To ensure that the sequence does not lie on a
|
|
Packit |
67cb25 |
subspace we force 32 of the entries to be linearly independent. We
|
|
Packit |
67cb25 |
take the 32 elements x[3], x[10], x[17], x[24], ..., 213 and apply
|
|
Packit |
67cb25 |
the following operations,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x[3] &= 11111111111111111111111111111111
|
|
Packit |
67cb25 |
x[3] |= 10000000000000000000000000000000
|
|
Packit |
67cb25 |
x[10] &= 01111111111111111111111111111111
|
|
Packit |
67cb25 |
x[10] |= 01000000000000000000000000000000
|
|
Packit |
67cb25 |
x[17] &= 00111111111111111111111111111111
|
|
Packit |
67cb25 |
x[17] |= 00100000000000000000000000000000
|
|
Packit |
67cb25 |
.... ...
|
|
Packit |
67cb25 |
x[206] &= 00000000000000000000000000000111
|
|
Packit |
67cb25 |
x[206] |= 00000000000000000000000000000100
|
|
Packit |
67cb25 |
x[213] &= 00000000000000000000000000000011
|
|
Packit |
67cb25 |
x[213] |= 00000000000000000000000000000010
|
|
Packit |
67cb25 |
x[220] &= 00000000000000000000000000000001
|
|
Packit |
67cb25 |
x[220] |= 00000000000000000000000000000001
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
i.e. if we consider the bits of the 32 elements as forming a 32x32
|
|
Packit |
67cb25 |
array then we are setting the diagonal bits of the array to one and
|
|
Packit |
67cb25 |
masking the lower triangle below the diagonal to zero.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
With this initialization procedure the theoretical value of
|
|
Packit |
67cb25 |
x_{10001} is 1100653588 for s = 1 (Actually I got this by running
|
|
Packit |
67cb25 |
the original code). The subscript 10001 means (1) seed the
|
|
Packit |
67cb25 |
generator with s = 1 and then do 10000 actual iterations.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The period of this generator is about 2^250.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The algorithm works for any number of bits. It is implemented here
|
|
Packit |
67cb25 |
for 32 bits.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
From: S. Kirkpatrick and E. Stoll, "A very fast shift-register
|
|
Packit |
67cb25 |
sequence random number generator", Journal of Computational Physics,
|
|
Packit |
67cb25 |
40, 517-526 (1981). */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline unsigned long int r250_get (void *vstate);
|
|
Packit |
67cb25 |
static double r250_get_double (void *vstate);
|
|
Packit |
67cb25 |
static void r250_set (void *state, unsigned long int s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
typedef struct
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int i;
|
|
Packit |
67cb25 |
unsigned long x[250];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
r250_state_t;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline unsigned long int
|
|
Packit |
67cb25 |
r250_get (void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
r250_state_t *state = (r250_state_t *) vstate;
|
|
Packit |
67cb25 |
unsigned long int k;
|
|
Packit |
67cb25 |
int j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int i = state->i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i >= 147)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
j = i - 147;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
j = i + 103;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
k = state->x[i] ^ state->x[j];
|
|
Packit |
67cb25 |
state->x[i] = k;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i >= 249)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
state->i = 0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
state->i = i + 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return k;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
r250_get_double (void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return r250_get (vstate) / 4294967296.0 ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
r250_set (void *vstate, unsigned long int s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
r250_state_t *state = (r250_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s == 0)
|
|
Packit |
67cb25 |
s = 1; /* default seed is 1 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->i = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define LCG(n) ((69069 * n) & 0xffffffffUL)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < 250; i++) /* Fill the buffer */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
s = LCG (s);
|
|
Packit |
67cb25 |
state->x[i] = s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Masks for turning on the diagonal bit and turning off the
|
|
Packit |
67cb25 |
leftmost bits */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
unsigned long int msb = 0x80000000UL;
|
|
Packit |
67cb25 |
unsigned long int mask = 0xffffffffUL;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < 32; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int k = 7 * i + 3; /* Select a word to operate on */
|
|
Packit |
67cb25 |
state->x[k] &= mask; /* Turn off bits left of the diagonal */
|
|
Packit |
67cb25 |
state->x[k] |= msb; /* Turn on the diagonal bit */
|
|
Packit |
67cb25 |
mask >>= 1;
|
|
Packit |
67cb25 |
msb >>= 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static const gsl_rng_type r250_type =
|
|
Packit |
67cb25 |
{"r250", /* name */
|
|
Packit |
67cb25 |
0xffffffffUL, /* RAND_MAX */
|
|
Packit |
67cb25 |
0, /* RAND_MIN */
|
|
Packit |
67cb25 |
sizeof (r250_state_t),
|
|
Packit |
67cb25 |
&r250_set,
|
|
Packit |
67cb25 |
&r250_get,
|
|
Packit |
67cb25 |
&r250_get_double};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const gsl_rng_type *gsl_rng_r250 = &r250_type;
|