|
Packit |
67cb25 |
/* rng/zuf.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 |
/* It is crucial that m == n-273 mod 607 at all times;
|
|
Packit |
67cb25 |
For speed of execution, however, this is never enforced.
|
|
Packit |
67cb25 |
Instead is is set in the initializer: note 607-273=334
|
|
Packit |
67cb25 |
Note also that the state.u[607] is not initialized */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline unsigned long int zuf_get (void *vstate);
|
|
Packit |
67cb25 |
static double zuf_get_double (void *vstate);
|
|
Packit |
67cb25 |
static void zuf_set (void *state, unsigned long int s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static const unsigned long int zuf_randmax = 16777216; /* 2^24 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
typedef struct
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int n;
|
|
Packit |
67cb25 |
unsigned long int u[607];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
zuf_state_t;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* The zufall package was implemented with float's, which is to say 24
|
|
Packit |
67cb25 |
bits of precision. Since I'm using long's instead, my RANDMAX
|
|
Packit |
67cb25 |
reflects this. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static inline unsigned long int
|
|
Packit |
67cb25 |
zuf_get (void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
zuf_state_t *state = (zuf_state_t *) vstate;
|
|
Packit |
67cb25 |
const int n = state->n;
|
|
Packit |
67cb25 |
const int m = (n - 273 + 607) % 607;
|
|
Packit |
67cb25 |
unsigned long int t = state->u[n] + state->u[m];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while (t > zuf_randmax)
|
|
Packit |
67cb25 |
t -= zuf_randmax;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->u[n] = t;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (n == 606)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
state->n = 0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
state->n = n + 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return t;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
zuf_get_double (void *vstate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return zuf_get (vstate) / 16777216.0 ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
zuf_set (void *vstate, unsigned long int s)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* A very elaborate seeding procedure is provided with the
|
|
Packit |
67cb25 |
zufall package; this is virtually a copy of that procedure */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Initialized data */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
long int kl = 9373;
|
|
Packit |
67cb25 |
long int ij = 1802;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Local variables */
|
|
Packit |
67cb25 |
long int i, j, k, l, m;
|
|
Packit |
67cb25 |
double x, y;
|
|
Packit |
67cb25 |
long int ii, jj;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
zuf_state_t *state = (zuf_state_t *) vstate;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
state->n = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* generates initial seed buffer by linear congruential */
|
|
Packit |
67cb25 |
/* method. Taken from Marsaglia, FSU report FSU-SCRI-87-50 */
|
|
Packit |
67cb25 |
/* variable seed should be 0 < seed <31328 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (s == 0)
|
|
Packit |
67cb25 |
s = 1802; /* default seed is 1802 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ij = s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
i = ij / 177 % 177 + 2;
|
|
Packit |
67cb25 |
j = ij % 177 + 2;
|
|
Packit |
67cb25 |
k = kl / 169 % 178 + 1;
|
|
Packit |
67cb25 |
l = kl % 169;
|
|
Packit |
67cb25 |
for (ii = 0; ii < 607; ++ii)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
x = 0.0;
|
|
Packit |
67cb25 |
y = 0.5;
|
|
Packit |
67cb25 |
/* 24 bits?? */
|
|
Packit |
67cb25 |
for (jj = 1; jj <= 24; ++jj)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
m = i * j % 179 * k % 179;
|
|
Packit |
67cb25 |
i = j;
|
|
Packit |
67cb25 |
j = k;
|
|
Packit |
67cb25 |
k = m;
|
|
Packit |
67cb25 |
l = (l * 53 + 1) % 169;
|
|
Packit |
67cb25 |
if (l * m % 64 >= 32)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
x += y;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
y *= 0.5;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
state->u[ii] = (unsigned long int) (x * zuf_randmax);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static const gsl_rng_type zuf_type =
|
|
Packit |
67cb25 |
{"zuf", /* name */
|
|
Packit |
67cb25 |
0x00ffffffUL, /* RAND_MAX */
|
|
Packit |
67cb25 |
0, /* RAND_MIN */
|
|
Packit |
67cb25 |
sizeof (zuf_state_t),
|
|
Packit |
67cb25 |
&zuf_set,
|
|
Packit |
67cb25 |
&zuf_get,
|
|
Packit |
67cb25 |
&zuf_get_double};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const gsl_rng_type *gsl_rng_zuf = &zuf_type;
|