|
Packit |
67cb25 |
/* rng/taus.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 maximally equidistributed combined Tausworthe
|
|
Packit |
67cb25 |
generator. The sequence is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x_n = (s1_n ^ s2_n ^ s3_n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s1_{n+1} = (((s1_n & 4294967294) <<12) ^ (((s1_n <<13) ^ s1_n) >>19))
|
|
Packit |
67cb25 |
s2_{n+1} = (((s2_n & 4294967288) << 4) ^ (((s2_n << 2) ^ s2_n) >>25))
|
|
Packit |
67cb25 |
s3_{n+1} = (((s3_n & 4294967280) <<17) ^ (((s3_n << 3) ^ s3_n) >>11))
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
computed modulo 2^32. In the three formulas above '^' means
|
|
Packit |
67cb25 |
exclusive-or (C-notation), not exponentiation. Note that the
|
|
Packit |
67cb25 |
algorithm relies on the properties of 32-bit unsigned integers (it
|
|
Packit |
67cb25 |
is formally defined on bit-vectors of length 32). I have added a
|
|
Packit |
67cb25 |
bitmask to make it work on 64 bit machines.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
We initialize the generator with s1_1 .. s3_1 = s_n MOD m, where
|
|
Packit |
67cb25 |
s_n = (69069 * s_{n-1}) mod 2^32, and s_0 = s is the user-supplied
|
|
Packit |
67cb25 |
seed.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The theoretical value of x_{10007} is 2733957125. The subscript
|
|
Packit |
67cb25 |
10007 means (1) seed the generator with s=1 (2) do six warm-up
|
|
Packit |
67cb25 |
iterations, (3) then do 10000 actual iterations.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The period of this generator is about 2^88.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
From: P. L'Ecuyer, "Maximally Equidistributed Combined Tausworthe
|
|
Packit |
67cb25 |
Generators", Mathematics of Computation, 65, 213 (1996), 203--213.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is available on the net from L'Ecuyer's home page,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme.ps
|
|
Packit |
67cb25 |
ftp://ftp.iro.umontreal.ca/pub/simulation/lecuyer/papers/tausme.ps
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Update: April 2002
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
There is an erratum in the paper "Tables of Maximally
|
|
Packit |
67cb25 |
Equidistributed Combined LFSR Generators", Mathematics of
|
|
Packit |
67cb25 |
Computation, 68, 225 (1999), 261--269:
|
|
Packit |
67cb25 |
http://www.iro.umontreal.ca/~lecuyer/myftp/papers/tausme2.ps
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
... the k_j most significant bits of z_j must be non-
|
|
Packit |
67cb25 |
zero, for each j. (Note: this restriction also applies to the
|
|
Packit |
67cb25 |
computer code given in [4], but was mistakenly not mentioned in
|
|
Packit |
67cb25 |
that paper.)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This affects the seeding procedure by imposing the requirement
|
|
Packit |
67cb25 |
s1 > 1, s2 > 7, s3 > 15.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The generator taus2 has been added to satisfy this requirement.
|
|
Packit |
67cb25 |
The original taus generator is unchanged.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Update: November 2002
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
There was a bug in the correction to the seeding procedure for s2.
|
|
Packit |
67cb25 |
It affected the following seeds 254679140 1264751179 1519430319
|
|
Packit |
67cb25 |
2274823218 2529502358 3284895257 3539574397 (s2 < 8).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline unsigned long int taus_get (void *vstate);
|
|
Packit |
67cb25 |
static double taus_get_double (void *vstate);
|
|
Packit |
67cb25 |
static void taus_set (void *state, unsigned long int s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
typedef struct
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
unsigned long int s1, s2, s3;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
taus_state_t;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline unsigned long
|
|
Packit |
67cb25 |
taus_get (void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
taus_state_t *state = (taus_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define MASK 0xffffffffUL
|
|
Packit |
67cb25 |
#define TAUSWORTHE(s,a,b,c,d) (((s &c) <<d) &MASK) ^ ((((s <<a) &MASK)^s) >>b)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->s1 = TAUSWORTHE (state->s1, 13, 19, 4294967294UL, 12);
|
|
Packit |
67cb25 |
state->s2 = TAUSWORTHE (state->s2, 2, 25, 4294967288UL, 4);
|
|
Packit |
67cb25 |
state->s3 = TAUSWORTHE (state->s3, 3, 11, 4294967280UL, 17);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return (state->s1 ^ state->s2 ^ state->s3);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
taus_get_double (void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return taus_get (vstate) / 4294967296.0 ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
taus_set (void *vstate, unsigned long int s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
taus_state_t *state = (taus_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s == 0)
|
|
Packit |
67cb25 |
s = 1; /* default seed is 1 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define LCG(n) ((69069 * n) & 0xffffffffUL)
|
|
Packit |
67cb25 |
state->s1 = LCG (s);
|
|
Packit |
67cb25 |
state->s2 = LCG (state->s1);
|
|
Packit |
67cb25 |
state->s3 = LCG (state->s2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* "warm it up" */
|
|
Packit |
67cb25 |
taus_get (state);
|
|
Packit |
67cb25 |
taus_get (state);
|
|
Packit |
67cb25 |
taus_get (state);
|
|
Packit |
67cb25 |
taus_get (state);
|
|
Packit |
67cb25 |
taus_get (state);
|
|
Packit |
67cb25 |
taus_get (state);
|
|
Packit |
67cb25 |
return;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
taus2_set (void *vstate, unsigned long int s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
taus_state_t *state = (taus_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s == 0)
|
|
Packit |
67cb25 |
s = 1; /* default seed is 1 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define LCG(n) ((69069 * n) & 0xffffffffUL)
|
|
Packit |
67cb25 |
state->s1 = LCG (s);
|
|
Packit |
67cb25 |
if (state->s1 < 2) state->s1 += 2UL;
|
|
Packit |
67cb25 |
state->s2 = LCG (state->s1);
|
|
Packit |
67cb25 |
if (state->s2 < 8) state->s2 += 8UL;
|
|
Packit |
67cb25 |
state->s3 = LCG (state->s2);
|
|
Packit |
67cb25 |
if (state->s3 < 16) state->s3 += 16UL;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* "warm it up" */
|
|
Packit |
67cb25 |
taus_get (state);
|
|
Packit |
67cb25 |
taus_get (state);
|
|
Packit |
67cb25 |
taus_get (state);
|
|
Packit |
67cb25 |
taus_get (state);
|
|
Packit |
67cb25 |
taus_get (state);
|
|
Packit |
67cb25 |
taus_get (state);
|
|
Packit |
67cb25 |
return;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static const gsl_rng_type taus_type =
|
|
Packit |
67cb25 |
{"taus", /* name */
|
|
Packit |
67cb25 |
0xffffffffUL, /* RAND_MAX */
|
|
Packit |
67cb25 |
0, /* RAND_MIN */
|
|
Packit |
67cb25 |
sizeof (taus_state_t),
|
|
Packit |
67cb25 |
&taus_set,
|
|
Packit |
67cb25 |
&taus_get,
|
|
Packit |
67cb25 |
&taus_get_double};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const gsl_rng_type *gsl_rng_taus = &taus_type;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static const gsl_rng_type taus2_type =
|
|
Packit |
67cb25 |
{"taus2", /* name */
|
|
Packit |
67cb25 |
0xffffffffUL, /* RAND_MAX */
|
|
Packit |
67cb25 |
0, /* RAND_MIN */
|
|
Packit |
67cb25 |
sizeof (taus_state_t),
|
|
Packit |
67cb25 |
&taus2_set,
|
|
Packit |
67cb25 |
&taus_get,
|
|
Packit |
67cb25 |
&taus_get_double};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const gsl_rng_type *gsl_rng_taus2 = &taus2_type;
|