From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 10015 invoked by alias); 7 Jan 2011 18:47:54 -0000 Received: (qmail 9971 invoked by uid 22791); 7 Jan 2011 18:47:53 -0000 X-SWARE-Spam-Status: No, hits=-1.0 required=5.0 tests=AWL,BAYES_50,T_RP_MATCHES_RCVD X-Spam-Check-By: sourceware.org Received: from mail.phy.duke.edu (HELO mail.phy.duke.edu) (152.3.182.2) by sourceware.org (qpsmtpd/0.43rc1) with SMTP; Fri, 07 Jan 2011 18:47:47 +0000 Received: from localhost (localhost [127.0.0.1]) by mail.phy.duke.edu (Postfix) with ESMTP id 9B79E780EA for ; Fri, 7 Jan 2011 13:45:49 -0500 (EST) Received: from mail.phy.duke.edu ([127.0.0.1]) by localhost (mail.phy.duke.edu [127.0.0.1]) (amavisd-new, port 10026) with LMTP id VSrC9+kp9UgR for ; Fri, 7 Jan 2011 13:45:49 -0500 (EST) Received: from lilith (unknown [10.180.136.65]) by mail.phy.duke.edu (Postfix) with ESMTP id 70822780DB for ; Fri, 7 Jan 2011 13:45:49 -0500 (EST) Date: Fri, 07 Jan 2011 18:47:00 -0000 From: "Robert G. Brown" To: GSL Discussion list Subject: Three random number generators... Message-ID: User-Agent: Alpine 2.02 (LFD 1266 2009-07-14) MIME-Version: 1.0 Content-Type: TEXT/PLAIN; format=flowed; charset=US-ASCII Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2011-q1/txt/msg00001.txt.bz2 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 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;