Three random number generators...

["Robert G. Brown" <rgb@phy.duke.edu>]

I have three random number generators that I think ought to go into the
GSL.  One is AES -- which is pretty self explanatory -- fairly fast,
cryptographic strength.

A second is Marsaglia's KISS (keep it simple, stupid) generator, which
is the fastest, best generator in dieharder at the moment -- I'm using a
GPL'd version of it by David Jones documented here:

   http://www.cs.ucl.ac.uk/staff/d.jones/GoodPracticeRNG.pdf

but using mt19937_1999 (in the GSL) to seed it.  It is possible to
improve it in a couple of ways indicated in the paper so that it is
"better" than a normalized uint as a float generator (in particular, so
it is truly random at bit granularity in floating point at the expense
of slowing it down, but I haven't tested how MUCH it will slow down yet.
It is also a good candidate for a 64 bit rng if/when the GSL is ready to
tackle 64 bit random numbers.  I append the code for it below, as it is
really easy and yet amazingly powerful (where it needs the usual lines
in the gsl rng includes).  Its only weakness is a relatively short period 
(compared to mt's) of 2^121 \approx 10^36.

However, this can easily be fixed by using a longer state array, and
Jones gives three versions of Marsaglia generators with absurdly long
periods (and excellent randomness).  I'm going to implement superkiss
next as it VECTORIZES the actual generation of the random numbers in a
way that might make it actually faster than any of the others in a
production environment (at the expense of an irrelevant slowdown if you
use <4096 rnds).

Finally, I've just implemented an XOR super-generator, that takes a list
(up to 100) of GSL or GSL-wrapped generators, seeds them from
mt19937_1999, and then xor's their output streams together, something
that can be no LESS random than the MOST random of the generators in the
stream, and something that makes the period (especially if any one of
the generators is mt19937) effectively infinite.  This obviously
(approximately linearly) slows as one adds more generators to the list,
but it gives me a very trustworthy "gold standard" rng to test dieharder
(where I don't care about speed, just positively identifying 
tests that fail when pushed too far).  It is GSL-wrapped, but I'm not
certain that the user interface is optimal for the GSL and it might be
better provided for people interested in it by looking at the dieharder
sources (which are not yet posted as I'm still cleaning it up and
testing it).

Basically, I'm suggesting that it is time to update the random number
generators in the GSL to make them more consistent with state of the art
in speed and tested randomness.

    rgb

Robert G. Brown	                       http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb@phy.duke.edu


/*
  *  rng_kiss.c
  *
  *  This is a version of (j)kiss (in C) from:
  *
  * Good Practice in (Pseudo) Random Number Generation for
  * Bioinformatics Applications
  *
  * http://www.cs.ucl.ac.uk/staff/d.jones/GoodPracticeRNG.pdf
  *
  * David Jones, UCL Bioinformatics Group
  * (E-mail: d.jones@cs.ucl.ac.uk)
  * (Last revised May 7th 2010)
  *
  *  GSL/GPL packaged for dieharder by Robert G. Brown 01/06/11
  *
  * David has kindly agreed to make this code vailable under the GPL so
  * it can be integrated with the dieharder package and/or the Gnu
  * Scentific Libraary (GSL).
  *
  */

#include <dieharder/libdieharder.h>

static unsigned long int kiss_get (void *vstate);
static double kiss_get_double (void *vstate);
static void kiss_set (void *vstate, unsigned long int s);

typedef struct {
  /*
   * Seed variables.  Note that kiss requires a moderately complex
   * seeding using a "seed rng" that we will arbitrarily set to be
   * mt19937_1999 from the GSL.
   */
  unsigned int x;
  unsigned int y;
  unsigned int z;
  unsigned int c;
} kiss_state_t;


static unsigned long int kiss_get (void *vstate)
{

  kiss_state_t *state = vstate;
  unsigned long long t;
  state->x = 314527869 * state->x + 1234567;
  state->y ^= state->y << 5;
  state->y ^= state->y >> 7;
  state->y ^= state->y << 22;
  t = 4294584393ULL * state->z + state->c;
  state->c = t >> 32;
  state->z = t;
  return (unsigned int)(state->x + state->y + state->z);

}

static double kiss_get_double (void *vstate)
{
   return kiss_get (vstate) / (double) UINT_MAX;
}

static void
kiss_set (void *vstate, unsigned long int s)
{

  kiss_state_t *state = (kiss_state_t *) vstate;

  uint seed_seed;
  gsl_rng *seed_rng;    /* random number generator used to seed kiss */

  /*
   * kiss needs four random number seeds.  They have to be reproducible
   * from a single seed in order to be consistent with the GSL.  We
   * therefore have to do a two step process where we use seed to
   * seed an existing GSL generator (mt19937_1999) and take the
   * first four returns as the seed for this generator.
   */
  seed_rng = gsl_rng_alloc(dh_rng_types[14]);
  seed_seed = s;
  gsl_rng_set(seed_rng,seed_seed);
  state->x = gsl_rng_get(seed_rng);
  while (!(state->y = gsl_rng_get(seed_rng))); /* y must not be zero! */
  state->z = gsl_rng_get(seed_rng);

  /*
   * We don't really need to set c as well but let's anyway.
   * Notes: offset c by 1 to avoid z=c=0; should be less than 698769069.
   */
  state->c = gsl_rng_get(seed_rng) % 698769068 + 1;

  return;

}

static const gsl_rng_type kiss_type =
{"kiss",			/* name */
  UINT_MAX,			/* RAND_MAX */
  0,				/* RAND_MIN */
  sizeof (kiss_state_t),
  &kiss_set,
  &kiss_get,
  &kiss_get_double};

const gsl_rng_type *gsl_rng_kiss = &kiss_type;