Suprangen: Sydney University Pseudorandom Number Generators


SourceForge.net project page
Suprangen project page hosted at: SourceForge.net Logo
About Suprangen
Suprangen: Sydney University Pseudorandom Number Generators This software is owned by the University of Sydney.

Release 0.4.2 2012-01-24

User Guide

1. Introduction

Suprangen (Sydney University Pseudo RAndom Number GENerators) is a library of pseudorandom number generators for the purpose of Monte Carlo simulation and similar applications.

It was initially created to provide a Fortran source code library of random number generators to enable Physics codes written in Fortran to be archived in such a way as to ensure repeatability, even in the face of changes to Fortran compilers. It has since grown in a number of ways to provide added value.

The main features of Suprangen are listed below and some of these features are covered in more detail in the following sections.

1.1. Simple static interface

For a simple program which uses one stream of uniform (0,1) random numbers, just call su_uniform_random(). For random numbers with standard normal distribution, call su_normal_random().

1.2. Multiple stream interface in C

The C interface uses pointers to structures corresponding to generators. These provide multiple independent streams of pseudo-random numbers. When used with a high-quality generator and careful seeding, the independent streams are practically uncorrelated.

1.3. Interfaces to GNU Scientific Library

The interfaces enable Suprangen to use any GSL generator, and enable GSL to use Suprangen generators. In particular, it is possible to use the GSL random distributions with Suprangen generators.

1.4. Fortran interface

The static interface, including su_uniform_random() and su_normal_random() is provided for Fortran 90 and above.

1.5. Fortran source code

Suprangen includes source code implementations of Xorgens XOR4096 [1] and Mersenne Twister MT19937 [3]. These produce identical output to the corresponding generators written in C.

1.6. Double precision uniform generators with 32 or 52 random bits

The uniform generator su_uniform_random() returns double precision results, based on a generator which supplies uniformly distributed 32 bit unsigned integers. Suprangen provides two different levels of precision, 32 bit, using a fast conversion algorithm inspired by Doornik [2], and 52 bit, using an algorithm inspired by [2] and Faustino de Sousa [3]. The 52 bit algorithm uses 5 calls to the 32 bit unsigned integer generator to return 3 double precision numbers.

1.7. Fast normal generator

The standard normal generator su_normal_random() uses the well-known, efficient modified polar form of the Box-Muller algorithm, [0,7] due to Marsaglia [5,6]. This algorithm generates two normal random numbers at one time. Unlike the GSL implementation, su_normal_random() uses the Polar Box-Muller algorithm only on every second call.

1.8. New, efficient, high quality 32 bit unsigned generators

As well as other 32 bit unsigned integer generators, Suprangen provides the 32 bit XOR4096 generator. XOR4096 is a a generator from the Xorgens collection [1] of improved xorshift generators.

1.9. Sensible defaults

The defaults are to use the Xorgens XOR4096 generator with 32 random bits of precision. XOR4096 is both fast and high quality.

1.10. Choice of algorithm delayed until linking

It is possible to write C or Fortran programs using a generic interface which is independent of which 32 bit unsigned generator is being used, and is independent of whether 32 or 52 random bits are generated. The specific choice of generator is delayed to linking time. This means you can have one set of source code and use different generators determined by the contents of your Makefile.

1.11. Extensively tested

The 32 and 52 bit uniform generators have been extensively tested using the test batteries from TestU01 [4]. In particular, both the 32 bit generator and the 52 bit generator using XOR4096 pass the stringent BigCrush test.

2. Suprangen types

2.1 Types of generator

The Suprangen library provides interfaces to three types of pseudorandom number generator. These are:

uniform uint32_t
A uniform uint32_t generator returns pseudorandom 32 bit unsigned integers uniformly distributed from 0 to 2^32-1.
uniform U(0,1)
A uniform U(0,1) generator returns pseudorandom double precision floating point numbers, uniformly distributed on the open interval (0,1).
normal N(0,1)
A normal N(0,1) generator returns pseudorandom double precision floating point numbers, normally distributed with mean 0 and variance 1.

These three types of generators are stacked, with the normal N(0,1) generator using a uniform U(0,1) generator, which in turn uses a uniform uint32_t generator.

3. Details of Suprangen features.

3.1. Simple static C interface

For a simple program which uses one stream of uniform U(0,1) random numbers, your source code needs to look like:

 #include <suprangen_static.h>
 /* ... */
 double x;
 /* ... */
 x = su_uniform_random();
 /* ... */

Other functions in the in the default static interface for the uniform generator are described in [M, pp. 62-64] in the section "8.11 include/su_static_uniform_gen.h File Reference", and in the file include/su_static_uniform_gen.h itself. Briefly, these are:

 void su_init_uniform_gen(uint32_t my_seed);
   /* Initializes the static default uniform generator from my_seed. */

 const char* su_uniform_name(void);
   /* Returns the name of the generator. */

 const char* su_uniform_version(void);
   /* Returns the version of the generator. */

 int su_uniform_bits(void);
   /* Returns the number of random bits returned by the generator. */

 uint32_t su_uniform_seed(void);
   /* Returns the current seed of the generator. */

For random numbers with standard normal distribution N(0,1), your source code needs to look like:

 #include <suprangen_static.h>
 /* ... */
 double x;
 /* ... */
 x = su_normal_random();
 /* ... */

Other functions in the in the default static interface for the normal N(0,1) generator are described in [M, pp. 59-61] in the section "8.9 include/su_static_normal_gen.h File Reference", and in the file include/su_static_normal_gen.h itself. Briefly, these are:

 void su_init_normal_gen(uint32_t seed);
   /* Initializes the static default normal generator from the given seed. */

 const char* su_normal_name(void);
   /* Returns the name of the generator. */

 const char* su_normal_version(void);
   /* Returns the version of the generator. */

 int su_normal_bits(void);
   /* Returns the number of random bits returned by the generator. */

 uint32_t su_normal_seed(void);
   /* Returns the current seed of the generator. */

Example

An example of the use of the simple static interface is given in the file test/test_suprangen.c. The corresponding Makefile is test/Makefile-c. With test as the current directory, the command:

 make -f Makefile-c
creates the file test_suprangen.exe. Running this file results in the output contained in test/test_suprangen.exe.out.

3.2. Multiple stream interface in C

The C interface uses pointers to structures corresponding to generators. These provide multiple independent streams of pseudo-random numbers. When used with a high-quality generator and careful seeding, the independent streams are practically uncorrelated.

To use the default multiple stream interface for uniform U(0,1) generators, your source code needs to look like:

 #include <su_uniform_gen.h>
 #include <su_uniform_default_gen.h>
 #include <su_uint32_default_gen.h>

 /* Create a default uniform generator using the given seed. */
 uint32_t my_seed = 0x12345678;
 uniform_gen_t* gen_p = su_new_uniform_gen(su_new_uint32_gen(my_seed));

 /* Return a uniform U(0,1) random double value. */
 double x = su_uniform_random_from(gen_p);

The idea is that gen_p represents a stream of uniform U(0,1) pseudorandom numbers seeded by my_seed. It is a pointer to a uniform_gen_t structure which is created when su_new_uniform_gen is called. The call to su_uniform_random_from obtains the next pseudorandom number from gen_p.

Other functions in the in the default multiple stream interface for uniform U(0,1) generators are described in [M, pp. 71-74] in the section "8.21 include/su_uniform_gen.h File Reference", and in the file include/su_uniform_gen.h itself. Briefly, these are:

 void su_set_uniform_gen(uniform_gen_t* gen_p, uint32_t seed);
   /* Seeds the generator. */

 const char* su_uniform_name_of(uniform_gen_t* gen_p);
   /* Returns the name of the generator. */

 const char* su_uniform_version_of(uniform_gen_t* gen_p);
   /* Returns the version. */

 int su_uniform_bits_of(uniform_gen_t* gen_p);
   /* Returns the number of random bits. */

 uint32_t su_uniform_seed_of(uniform_gen_t* gen_p);
   /* Returns the current seed. */

 uniform_gen_t* su_copy_uniform_gen(uniform_gen_t* gen_p);
   /* Copies state information from the generator. */

 void su_end_uniform_gen(uniform_gen_t* gen_p);
   /* Deletes state information from the generator. */

To use the default multiple stream interface for standard normal N(0,1) generators, your source code needs to look like:

 #include <su_normal_gen.h>
 #include <su_normal_default_gen.h>
 #include <su_uniform_default_gen.h>
 #include <su_uint32_default_gen.h>

 /* Create a default normal generator using the given seed. */
 uint32_t my_seed = 0x12345678;
 normal_gen_t*
   gen_p = su_new_normal_gen(su_new_uniform_gen(su_new_uint32_gen(my_seed)));

 /* Return a normal N(0,1) random double value. */
 double x = su_normal_random_from(gen_p);

The idea is that gen_p represents a stream of normal N(0,1) pseudorandom numbers seeded by my_seed. It is a pointer to a normal_gen_t structure which is created when su_new_normal_gen is called. The call to su_normal_random_from obtains the next pseudorandom number from gen_p.

Other functions in the in the default multiple stream interface for normal N(0,1) generators are described in [M, pp. 57-59] in the section "8.8 include/su_normal_gen.h File Reference", and in the file include/su_uniform_gen.h itself. Briefly, these are:

 void su_set_normal_gen(normal_gen_t* gen_p, uint32_t seed);
   /* Seeds the generator. */

 const char* su_normal_name_of(normal_gen_t* gen_p);
   /* Returns the name of the generator. */

 const char* su_normal_version_of(normal_gen_t* gen_p);
   /* Returns the version. */

 int su_normal_bits_of(normal_gen_t* gen_p);
   /* Returns the number of random bits. */

 uint32_t su_normal_seed_of(normal_gen_t* gen_p);
   /* Returns the current seed. */

 normal_gen_t* su_copy_normal_gen(normal_gen_t* gen_p);
   /* Copies state information from the generator. */

 void su_end_normal_gen(normal_gen_t* gen_p);
   /* Deletes state information from the generator. */

3.3. Interfaces to GNU Scientific Library

The interfaces enable Suprangen to use any GSL generator, and enable GSL to use Suprangen generators.

For documentation on the GSL RNG interface, see http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Generation.html

3.3.1 Suprangen using GSL generators

The interface for Suprangen to use GSL generators is provided via uint32_t generators. There are a number of ways to create a Suprangen uint32_t generator from a GSL generator:

1. From an existing GSL generator:

 #include <su_uint32_gsl_gen.h>
 uint32_gen_t* uint32_gen_p;

 /* Create a GSL generator using gsl_rng_ranlxd2 */
 const gsl_rng* gsl_rng_p = gsl_rng_alloc(gsl_rng_ranlxd2);
 /*...*/
 /* Create a new uint32_t generator structure for the existing GSL generator,
    seeded by my_seed.  */
 uint32_t my_seed = 0x12345678;
 uint32_gen_p = su_new_uint32_gsl_gen_from_gsl_rng(gsl_rng_p, my_seed);

2. From a GSL generator type:

 #include <su_uint32_gsl_gen.h>
 uint32_gen_t* uint32_gen_p;

 /* Create a new uint32_t generator structure using gsl_rng_ranlxd2,
    seeded by my_seed.  */
 uint32_t my_seed = 0x12345678;
 uint32_gen_p = su_new_uint32_gsl_gen(gsl_rng_ranlxd2, my_seed);

3. From the environment variable GSL_RNG_TYPE

In csh:

 setenv GSL_RNG_TYPE ranlxd2
In your program:
 #include <su_uint32_gsl_gen.h>
 uint32_gen_t* uint32_gen_p;

 /* Create a new uint32_t generator structure using gsl_rng_default,
    seeded by my_seed.  */
 uint32_t my_seed = 0x12345678;
 gsl_rng_env_setup();
 uint32_gen_p = su_new_uint32_gsl_gen(gsl_rng_default, my_seed);

3.3.2. GSL using Suprangen generators

In particular, it is possible to use the GSL random distributions with Suprangen generators.

3.4. Fortran interface

The static interface, including su_uniform_random() and su_normal_random() is provided for Fortran 90 and above.

3.5. Fortran source code

Suprangen includes source code implementations of Xorgens XOR4096 [1] and Mersenne Twister MT19937 [3]. These produce identical output to the corresponding generators written in C.

3.6. Double precision uniform generators with 32 or 52 random bits

The uniform generator su_uniform_random() returns double precision results, based on a generator which supplies uniformly distributed 32 bit unsigned integers. Suprangen provides two different levels of precision, 32 bit, using a fast conversion algorithm inspired by Doornik [2], and 52 bit, using an algorithm inspired by [2] and Faustino de Sousa [3]. The 52 bit algorithm uses 5 calls to the 32 bit unsigned integer generator to return 3 double precision numbers.

3.7. Fast normal generator

The standard normal generator su_normal_random() uses the well-known, efficient modified polar form of the Box-Muller algorithm, [0,7] due to Marsaglia [5,6]. This algorithm generates two normal random numbers at one time. Unlike the GSL implementation, su_normal_random() uses the Polar Box-Muller algorithm only on every second call.

3.8. New, efficient, high quality 32 bit unsigned generators

As well as other 32 bit unsigned integer generators, Suprangen provides the 32 bit XOR4096 generator. XOR4096 is a a generator from the Xorgens collection [1] of improved xorshift generators.

3.9. Sensible defaults

The defaults are to use the Xorgens XOR4096 generator with 32 random bits of precision. XOR4096 is both fast and high quality.

3.10. Choice of algorithm delayed until linking

It is possible to write C or Fortran programs using a generic interface which is independent of which 32 bit unsigned generator is being used, and is independent of whether 32 or 52 random bits are generated. The specific choice of generator is delayed to linking time. This means you can have one set of source code and use different generators determined by the contents of your Makefile.

3.11. Extensively tested

The 32 and 52 bit uniform generators have been extensively tested using the test batteries from TestU01 [4]. In particular, both the 32 bit generator and the 52 bit generator using XOR4096 pass the stringent BigCrush test.

References

[M] Suprangen Reference Manual, Version 0.4.2, 2012-01-24

[0] G. E. P. Box, M. E. Muller, "A note on the generation of random normal deviates", Annals Math. Stat, V. 29, No. 2, 1958, pp. 610-611.

[1] R. P. Brent, "Some long-period random number generators using shifts and xors", ANZIAM Journal 48 (CTAC2006)(1), C188?C202 (2007). Presented at the 13th Biennial Computational Techniques and Applications Conference (CTAC06), Townsville, 2-5 July 2006. arXiv:1004.3115v1 http://maths.anu.edu.au/~brent/random.html

[2] J. A. Doornik, "Conversion of High-Period Random Numbers to Floating Point".

[3] J. R. Faustino de Sousa, Module mt19937: Code converted to Fortran 95 by Jose Rui Faustino de Sousa, 2002-02-01. Enhanced version of this code is available at http://homepage.esoterica.pt/~jrfsousa/files/mt95.htm Based on [8].

[4] P. L'Ecuyer, R. Simard "Empirical Testing of Random Number Generators".

[5] G. Marsaglia, "Improving the polar method for generating a pair of random variables", Boeing Sci. res. Lab. D1-82-0203, 1962.

[6] G. Marsaglia, T.A. Bray, "A convenient method for generating normal variables", SIAM Review, V. 6, No. 2, July 1964, pp. 260-264.

[7] M.E. Muller, "A comparison of methods for generating normal deviates on digital computers", Journal of the ACM, V. 6, 1959, pp. 376-383.

[8] T. Nishimura, M. Matsumoto, "A C-program for MT19937, with initialization improved" 2002-01-26.


Presentations


Publications and preprints
Home page



Page Generated: Tuesday 24 January 2012