Monte Carlo in GSL

Robert G. Brown rgb@phy.duke.edu
Wed Feb 18 22:32:00 GMT 2004


On Wed, 18 Feb 2004, Przemyslaw Sliwa wrote:

> The question is: Do I have to execute the statement
> 
>   gsl_rng_env_setup();
> 
>   T = gsl_rng_default;
>   r = gsl_rng_alloc (T);
> 
> each time I compute the value of my function or just once at the beginning
> of the process where I solve the equation.

This stuff can go just once at the beginning of the program.

> My problem is: I want that the random numbers I generate are purelly
> random and not that the series of the numbers repeats.

"Purely random" is a (well-known) oxymoron as far as uniform deviates
are concerned, but "pseudo random" with good properites of randomness,
that might be doable:-)

If you want the rands to be "random" as a single sequence this is all
you need to do.  If you want a different sequence each time you run the
code, you need to seed the generator with a different seed each time you
run the code.  Picking a "good" seed isn't terribly easy, since at least
some rng's (admittedly probably not the default GSL generator) show
significant sequence-to-sequence correlations when started with nearby
seeds.  Seeding from e.g. the clock is thus a deprecated technique.

The "best" way to seed if you want all runs to be unique and for all the
runs to be maximally decorrelated at the seed level is to seed from the
kernel-based entropy generator in /dev/random (if it exists on your
system).  Save each seed as you use it, or (better yet) save your
results using the seed as the unique part of the label.  That way if by
any chance you get the same seed more than once (which will generate
identical runs) each duplicate run will just overwrite the results of
the previous run, and counting results files gives you an accurate count
of the number of unique, identically distributed runs.  Note that this
probability is VERY low and you can PROBABLY ignore it unless you plan
to make a LOT of different runs.

The following routine will generate a seed with /dev/random (if
possible).  You can arrange to fall back on using e.g. the clock if it
isn't.  Beware doing this on a cluster, where lots of jobs might start
on the same clock tick!

Then you can reseed with e.g.

  unsigned long int seed;
  seed = random_seed(); 
  gsl_rng_set(r,seed);

You can reseed AT ANY TIME to alter the pseudorandom sequence you are
using.  Be aware that reading /dev/random is SLOW and will HANG THE
PROCESS if there is insufficient entropy in the system to generate the
next number, until there is.  So reseed sparingly.  The default GSL
random number generator is "quite good" (fast and highly uncorrelated by
many measures) so you SHOULD need to seed just once...

HTH

   rgb

unsigned long int random_seed()
{

 unsigned int seed;
 FILE *devrandom;

 if ((devrandom = fopen("/dev/random","r")) == NULL) {
   fprintf(stderr,"Cannot open /dev/random, setting seed to 0\n");
   seed = 0;
 } else {
   fread(&seed,sizeof(seed),1,devrandom);
   if(verbose == D_SEED) printf("Got seed %u from /dev/random\n",seed);
   fclose(devrandom);
 }

 return(seed);

}



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





More information about the Gsl-discuss mailing list