|
Packit |
67cb25 |
/* rng/knuthran2002.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2007 Brian Gough
|
|
Packit |
67cb25 |
* Copyright (C) 2001, 2007 Brian Gough, Carlo Perassi
|
|
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 |
/*
|
|
Packit |
67cb25 |
* This generator is taken from
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Donald E. Knuth, The Art of Computer Programming, Volume 2, Section 3.6
|
|
Packit |
67cb25 |
* Third Edition, Addison-Wesley,
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* The modifications introduced in the 9th printing (2002) are
|
|
Packit |
67cb25 |
* included here; there's no backwards compatibility with the
|
|
Packit |
67cb25 |
* original. [ see http://www-cs-faculty.stanford.edu/~knuth/taocp.html ]
|
|
Packit |
67cb25 |
*
|
|
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 |
#define BUFLEN 1009 /* length of the buffer aa[] */
|
|
Packit |
67cb25 |
#define KK 100 /* the long lag */
|
|
Packit |
67cb25 |
#define LL 37 /* the short lag */
|
|
Packit |
67cb25 |
#define MM (1L << 30) /* the modulus */
|
|
Packit |
67cb25 |
#define TT 70 /* guaranteed separation between streams */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define is_odd(x) ((x) & 1) /* the units bit of x */
|
|
Packit |
67cb25 |
#define mod_diff(x, y) (((x) - (y)) & (MM - 1)) /* (x - y) mod MM */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline void ran_array (long int aa[], unsigned int n,
|
|
Packit |
67cb25 |
long int ran_x[]);
|
|
Packit |
67cb25 |
static inline unsigned long int ran_get (void *vstate);
|
|
Packit |
67cb25 |
static double ran_get_double (void *vstate);
|
|
Packit |
67cb25 |
static void ran_set (void *state, unsigned long int s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
typedef struct
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
unsigned int i;
|
|
Packit |
67cb25 |
long int aa[BUFLEN];
|
|
Packit |
67cb25 |
long int ran_x[KK]; /* the generator state */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
ran_state_t;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline void
|
|
Packit |
67cb25 |
ran_array (long int aa[], unsigned int n, long int ran_x[])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
unsigned int i;
|
|
Packit |
67cb25 |
unsigned int j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < KK; j++)
|
|
Packit |
67cb25 |
aa[j] = ran_x[j];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (; j < n; j++)
|
|
Packit |
67cb25 |
aa[j] = mod_diff (aa[j - KK], aa[j - LL]);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < LL; i++, j++)
|
|
Packit |
67cb25 |
ran_x[i] = mod_diff (aa[j - KK], aa[j - LL]);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (; i < KK; i++, j++)
|
|
Packit |
67cb25 |
ran_x[i] = mod_diff (aa[j - KK], ran_x[i - LL]);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline unsigned long int
|
|
Packit |
67cb25 |
ran_get (void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ran_state_t *state = (ran_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
unsigned int i = state->i;
|
|
Packit |
67cb25 |
unsigned long int v;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* fill buffer with new random numbers */
|
|
Packit |
67cb25 |
ran_array (state->aa, BUFLEN, state->ran_x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
v = state->aa[i];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->i = (i + 1) % KK;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return v;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
ran_get_double (void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ran_state_t *state = (ran_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return ran_get (state) / 1073741824.0; /* RAND_MAX + 1 */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
ran_set (void *vstate, unsigned long int s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ran_state_t *state = (ran_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
long x[KK + KK - 1]; /* the preparation buffer */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
register int j;
|
|
Packit |
67cb25 |
register int t;
|
|
Packit |
67cb25 |
register long ss;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s == 0 )
|
|
Packit |
67cb25 |
s = 314159; /* default seed used by Knuth */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ss = (s + 2)&(MM-2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < KK; j++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
x[j] = ss; /* bootstrap the buffer */
|
|
Packit |
67cb25 |
ss <<= 1;
|
|
Packit |
67cb25 |
if (ss >= MM) /* cyclic shift 29 bits */
|
|
Packit |
67cb25 |
ss -= MM - 2;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
x[1]++; /* make x[1] (and only x[1]) odd */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ss = s & (MM - 1);
|
|
Packit |
67cb25 |
t = TT - 1;
|
|
Packit |
67cb25 |
while (t)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = KK - 1; j > 0; j--) /* square */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
x[j + j] = x[j];
|
|
Packit |
67cb25 |
x[j + j - 1] = 0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = KK + KK - 2; j >= KK; j--)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
x[j - (KK - LL)] = mod_diff (x[j - (KK - LL)], x[j]);
|
|
Packit |
67cb25 |
x[j - KK] = mod_diff (x[j - KK], x[j]);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (is_odd (ss))
|
|
Packit |
67cb25 |
{ /* multiply by "z" */
|
|
Packit |
67cb25 |
for (j = KK; j > 0; j--)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
x[j] = x[j - 1];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
x[0] = x[KK]; /* shift the buffer cyclically */
|
|
Packit |
67cb25 |
x[LL] = mod_diff (x[LL], x[KK]);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (ss)
|
|
Packit |
67cb25 |
ss >>= 1;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
t--;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < LL; j++)
|
|
Packit |
67cb25 |
state->ran_x[j + KK - LL] = x[j];
|
|
Packit |
67cb25 |
for (; j < KK; j++)
|
|
Packit |
67cb25 |
state->ran_x[j - LL] = x[j];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j< 10; j++)
|
|
Packit |
67cb25 |
ran_array(x, KK+KK-1, state->ran_x); /* warm things up */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->i = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static const gsl_rng_type ran_type = {
|
|
Packit |
67cb25 |
"knuthran2002", /* name */
|
|
Packit |
67cb25 |
0x3fffffffUL, /* RAND_MAX = (2 ^ 30) - 1 */
|
|
Packit |
67cb25 |
0, /* RAND_MIN */
|
|
Packit |
67cb25 |
sizeof (ran_state_t),
|
|
Packit |
67cb25 |
&ran_set,
|
|
Packit |
67cb25 |
&ran_get,
|
|
Packit |
67cb25 |
&ran_get_double
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const gsl_rng_type *gsl_rng_knuthran2002 = &ran_type;
|