/* rand31-park-miller-carta-int.c Version 1.00 2005 September 21
*
* Robin Whittle rw@firstpr.com.au
*
* Double-precision floating point implementation of the Park Miller
* "minimal standard" linear congruential pseudo-random number
* generator.
*
* This file and its .h file is intended to be used in other projects.
*
* The accompanying files rand31-park-miller-carta-int.c/h have a fast
* implementation of the Park Miller (1988) linear congruential
* pseudo-random number generator, using David G. Carta's optimisation
* which needs only 32 bit integer math, and no division.
*
* A test program enables the speed of each approach to be tested by
* making each one run through the entire pseudo-random sequence once:
*
* rand31-park-miller-carta-c-test.c
*
* On an 800MHz Pentium III, with GCC and no optimisations, the integer
* version produced 13 million results a second, running in a simple
* test loop, and the floating point version produced 3.6 million.
*
* C++ versions are also available.
*
* References:
*
* Stephen K. Park and Keith W. Miller
* Random Number Generators: Good Ones are Hard to Find
* Communications of the ACM, Oct 1988, Vol 31 Number 10 1192-1201
*
* Like the other two papers, this one is normally only available
* from the ACM site via subscription. You should be able to
* access this paper electronically or in print at a university
* library. You may also be able to find the .PDF wild on the
* Net. Search for "p1192-park.pdf". For instance:
*
* http://www-scf.usc.edu/~csci105/links/p1192-park.pdf
*
* David F. Carta
* Two Fast Implementations of the "Minimal Standard" Random Number Generator
* Communications of the ACM, Jan 1990, Vol 33 Number 1 87-88 (p87-carta.pdf)
*
* George Marsaglia; Stephen J. Sullivan; Stephen K. Park, Keith W. Miller,
* Paul K. Stockmeyer
* Remarks on Choosing and Implementing Random Number Generators
* Communications of the ACM, Jul 1993, Vol 36 Number 7 105-110 (p105-crawford.pdf)
*
* The following code is public domain. If you use this code, I request that
* you keep the comments with it, to save some poor soul from having to figure
* out the history behind it. If you use a PRNG, you should research its
* pedigree.
*
* Copyright public domain Robin Whittle 2005
*
* For a full explanation, latest updates and the history of these
* algorithms, see:
*
* http://www.firstpr.com.au/dsp/rand31/
*
* When compiling into the test program I use:
*
* gcc rand31-park-miller-carta-c-test.c -o rand31-pmc-c-test -lm
*
* The -lm was necessary to stop the compiler complaining about fmod.
*
*/
#include "rand31-park-miller-fp.h"
#include
/* rand31pm_next()
*
* 31 bit Pseudo Random Number Generator
* based on Park Miller "Integer
* Version 1" - but done with double-
* precision floating point so we are not
* concerned with the limits of integer
* operations. This is not intended for
* fast operation - but maybe it is
* faster than the integer implementation
* on some CPUs.
*
* Generate next pseudo random number.
*
* Multiplier constant = 16807 = 7^5.
* This is 15 bits.
*
* Park and Miller in 1993 recommend
* 48271, which they say produces a
* somewhat better quality of
* pseudo-random results.
*/
#define constapm 16807
/* #define constapm 48271
*/
/* Modulus constant = 2^31 - 1 =
* 0x7FFFFFFFF. Use .0 to deter compiler
* from complaining about a very large
* integer constant.
*/
#define constmpm 2147483647.0
long unsigned int rand31pm_next()
{
double const a = constapm;
double const m = constmpm;
/* This is the linear congrentual
* generator:
*
* Multiply the old seed by constant a
* and take the modulus of the result
* (the remainder of a division) by
* constant m.
*/
return (seed31pm = (long)(fmod((seed31pm * a), m)) );
}
/*---------------------------------------------------------------------------*/
/* rand31pm_seedi()
*
* Set the seed from a long unsigned
* integer. If zero is used, then
* the seed will be set to 1.
*/
void rand31pm_seedi(long unsigned int seedin)
{
if (seedin == 0) seedin = 1;
seed31pm = seedin;
}
/* rand31pm_ranlui()
*
* Return next pseudo-random value as
* a long unsigned integer.
*/
long unsigned int rand31pm_ranlui(void)
{
return rand31pm_next();
}
/* rand31pm_ranf()
*
* Return next pseudo-random value as
* a floating point value.
*/
float rand31pm_ranf(void)
{
return (rand31pm_next() / 2147483647.0 );
}