diff options
Diffstat (limited to 'lib/util/rand.c')
| -rw-r--r-- | lib/util/rand.c | 294 |
1 files changed, 294 insertions, 0 deletions
diff --git a/lib/util/rand.c b/lib/util/rand.c new file mode 100644 index 00000000..6a0a2cdb --- /dev/null +++ b/lib/util/rand.c @@ -0,0 +1,294 @@ +/* + +Pseudo-random number generator for Netpbm + +The interface provided herein should be flexible enough for anybody +who wishes to use some other random number generator. + +--- + +If you desire to implement a different generator, or writing an original +one, first take a look at the random number generator section of the +GNU Scientific Library package (GSL). + +GNU Scientific Library +https://www.gnu.org/software/gsl/ + +GSL Random Number Generators +https://wnww.gnu.org/software/gsl/doc/html/rng.html + +*/ + +#include <stdio.h> +#include <stdlib.h> +#include <inttypes.h> +#include <strings.h> +#include <time.h> +#include <float.h> +#include <math.h> + +#include "netpbm/pm_c_util.h" +#include "netpbm/mallocvar.h" +#include "netpbm/pm.h" +#include "netpbm/rand.h" + +/*----------------------------------------------------------------------------- + Use +------------------------------------------------------------------------------- + Typical usage: + + #include "rand.h" + + ... + + myfunction( ... , unsigned int const seed , ... ) { + + struct randSt; + + ... + + pm_randinit(&randSt); + pm_srand(&randSt, seed); // pm_srand2() is often more useful + + ... + + pm_rand(&randSt); + + ... + + pm_randterm(&randSt); + + } +-----------------------------------------------------------------------------*/ + + + +/*----------------------------------------------------------------------------- + Design note +------------------------------------------------------------------------------- + + Netpbm code contains multiple random number generators. Stock Netpbm always + uses an internal pseudo-random number generator that implements the Mersenne + Twister method and does not rely on any randomness facility of the operating + system, but it is easy to compile an alternative version that uses others. + + The Mersenne Twister method was new to Netpbm in Netpbm 10.94 (March 2021). + Before that, Netpbm used standard OS-provided facilities. + + Programs that use random numbers have existed in Netpbm since PBMPlus days. + The system rand() function was used in instances randomness was required; + exceptions were rare and all of them appear to be errors on the part of the + original author. + + Although the rand() function is available in every system on which Netpbm + runs, differences exist in the underlying algorithm, so that Netpbm programs + produce different output on different systems even when the user specifies + the same random number seed. + + This was not considered a problem in the early days. Deterministic + operation was not a feature users requested and it was impossible regardless + of the random number generation method on most programs because they did not + allow a user to specify a seed for the generator. + + This state of affairs changed as Netpbm got firmly established as a + base-level system package. Security became critical for many users. A + crucial component of quality control is automated regression tests (="make + check"). Unpredictable behavior gets in the way of testing. One by one + programs were given the -randomseed (or -seed) option to ensure reproducible + results. Often this was done as new test cases were written. However, + inconsistent output caused by system-level differences in rand() + implementation remained a major obstacle. + + In 2020 the decision was made to replace all calls to rand() in the Netpbm + source code with an internal random number generator. We decided to use the + Mersenne Twister, which is concise, enjoys a fine reputation and is + available under liberal conditions (see below.) +-----------------------------------------------------------------------------*/ + + +void +pm_srand(struct pm_randSt * const randStP, + unsigned int const seed) { +/*---------------------------------------------------------------------------- + Initialize (or "seed") the random number generation sequence with value + 'seed'. +-----------------------------------------------------------------------------*/ + pm_randinit(randStP); + + randStP->vtable.srand(randStP, seed); + + randStP->seed = seed; +} + + + +void +pm_srand2(struct pm_randSt * const randStP, + bool const seedValid, + unsigned int const seed) { +/*---------------------------------------------------------------------------- + Seed the random number generator. If 'seedValid' is true, use 'seed".. + Otherwise, use pm_randseed(). + + For historical reasons pm_randseed() is defined in libpm.c rather than + this source file. +-----------------------------------------------------------------------------*/ + pm_srand(randStP, seedValid ? seed : pm_randseed() ); + +} + + + +unsigned long int +pm_rand(struct pm_randSt * const randStP) { +/*---------------------------------------------------------------------------- + An integer random number in the interval [0, randStP->max]. +-----------------------------------------------------------------------------*/ + return randStP->vtable.rand(randStP); +} + + + +double +pm_drand(struct pm_randSt * const randStP) { +/*---------------------------------------------------------------------------- + A floating point random number in the interval [0, 1). + + Although the return value is declared as double, the actual value will have + no more precision than a single call to pm_rand() provides. This is 32 bits + for Mersenne Twister. +-----------------------------------------------------------------------------*/ + return (double) pm_rand(randStP) / randStP->max; +} + + + +void +pm_gaussrand2(struct pm_randSt * const randStP, + double * const r1P, + double * const r2P) { +/*---------------------------------------------------------------------------- + Generate two Gaussian (or normally) distributed random numbers *r1P and + *r2P. + + Mean = 0, Standard deviation = 1. + + This is called the Box-Muller method. + + For details of this algorithm and other methods for producing + Gaussian random numbers see: + + http://www.doc.ic.ac.uk/~wl/papers/07/csur07dt.pdf +-----------------------------------------------------------------------------*/ + double u1, u2; + + u1 = pm_drand(randStP); + u2 = pm_drand(randStP); + + if (u1 < DBL_EPSILON) + u1 = DBL_EPSILON; + + *r1P = sqrt(-2.0 * log(u1)) * cos(2.0 * M_PI * u2); + *r2P = sqrt(-2.0 * log(u1)) * sin(2.0 * M_PI * u2); +} + + + +double +pm_gaussrand(struct pm_randSt * const randStP) { +/*---------------------------------------------------------------------------- + A Gaussian (or normally) distributed random number. + + Mean = 0, Standard deviation = 1. + + If a randStP->gaussCache has a value, return that value. Otherwise call + pm_gaussrand2; return one generated value, remember the other. +-----------------------------------------------------------------------------*/ + double retval; + + if (!randStP->gaussCacheValid) { + pm_gaussrand2(randStP, &retval, &randStP->gaussCache); + randStP->gaussCacheValid = true; + } else { + retval = randStP->gaussCache; + randStP->gaussCacheValid = false; + } + + return retval; +} + + + +uint32_t +pm_rand32(struct pm_randSt * const randStP) { +/*----------------------------------------------------------------------------- + Generate a 32-bit random number. + + This is a provision for users who select a non-default random number + generator which returns less than 32 bits per call. Many system generators + are known to return 31 bits (max = 2147483647 or 0x7FFFFFFF). + + This does not work with generators that return less than 11 bits per call. + The least we know of is the archaic RANDU, which generates 15 bits (max = + 32767 or 0x7FFF). +-----------------------------------------------------------------------------*/ + unsigned int const randMax = randStP->max; + + uint32_t retval; + + if (randMax >= 0xFFFFFFFF) + retval = pm_rand(randStP); + else { + uint32_t scale; + + retval = 0; /* Initial value */ + + for (scale = 0xFFFFFFFF; scale > 0; scale /= (randMax +1)) + retval *= (randMax + 1) + pm_rand(randStP); + } + + return retval;; +} + + + +void +pm_randinit(struct pm_randSt * const randStP) { +/*---------------------------------------------------------------------------- + Initialize the random number generator. +-----------------------------------------------------------------------------*/ + switch (PM_RANDOM_NUMBER_GENERATOR) { + case PM_RAND_SYS_RAND: + randStP->vtable = pm_randsysrand_vtable; + break; + case PM_RAND_SYS_RANDOM: + randStP->vtable = pm_randsysrandom_vtable; + break; + case PM_RAND_MERSENNETWISTER: + randStP->vtable = pm_randmersenne_vtable; + break; + default: + pm_error("INTERNAL ERROR: Invalid value of " + "PM_RANDOM_NUMBER_GENERATOR (random number generator " + "engine type): %u", PM_RANDOM_NUMBER_GENERATOR); + } + + randStP->vtable.init(randStP); + + randStP->gaussCacheValid = false; +} + + + +void +pm_randterm(struct pm_randSt * const randStP) { +/*---------------------------------------------------------------------------- + Tear down the random number generator. +-----------------------------------------------------------------------------*/ + if (randStP->stateP) + free(randStP->stateP); +} + + + + |