This is the mail archive of the
`gsl-discuss@sourceware.org`
mailing list for the GSL project.

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |

Other format: | [Raw text] |

*From*: "Robert G. Brown" <rgb at phy dot duke dot edu>*To*: GSL Discussion list <gsl-discuss at sources dot redhat dot com>*Date*: Fri, 7 Jan 2011 13:47:41 -0500 (EST)*Subject*: Three random number generators...

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:

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).

(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.

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). * */

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) {

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;

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};

**Follow-Ups**:**Re: Three random number generators...***From:*Brian Gough

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |