diff options
Diffstat (limited to 'lib/util/rand.c')
| -rw-r--r-- | lib/util/rand.c | 320 |
1 files changed, 320 insertions, 0 deletions
diff --git a/lib/util/rand.c b/lib/util/rand.c new file mode 100644 index 00000000..fd083c3b --- /dev/null +++ b/lib/util/rand.c @@ -0,0 +1,320 @@ +/* + +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 instances where random numbers are used. + 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 the operating system function, or some + other generator. + + 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 where 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 + produced different output on different systems even when the user specified + 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; +} + + + +double +pm_drand1(struct pm_randSt * const randStP) { +/*---------------------------------------------------------------------------- + A floating point random number in the interval [0, 1). +-----------------------------------------------------------------------------*/ + return (double) pm_rand(randStP) / ((double) randStP->max + 1.0); +} + + + +double +pm_drand2(struct pm_randSt * const randStP) { +/*---------------------------------------------------------------------------- + A floating point random number in the interval (0, 1). +-----------------------------------------------------------------------------*/ + return ((double) pm_rand(randStP) + 0.5) / ((double) randStP->max + 1.0); +} + + + +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_drand2(randStP); + u2 = pm_drand1(randStP); + + *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. +-----------------------------------------------------------------------------*/ + unsigned int const randMax = randStP->max; + + /* 'randMax is a power of 2 minus 1. Function pm_randinit() rejects + generators which do not satisfy this condition. It is unlikely that + such odd generators actually exist. + + Many system generators are known to return 31 bits (max = 2147483647 or + 0x7FFFFFFF). Historically, there were generators that returned only 15 + bits. + */ + + 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 = 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); + + if (randStP->max == 0) + pm_error("Random number generator maximum value must be positive."); + else if (((long int) randStP->max & (long int) (randStP->max + 1)) != 0x0L) + pm_error("Non-standard random number generator with maximum value " + "%u. Cannot handle maximum values which are not powers " + "of 2 minus 1", randStP->max); + + randStP->gaussCacheValid = false; +} + + + +void +pm_randterm(struct pm_randSt * const randStP) { +/*---------------------------------------------------------------------------- + Tear down the random number generator. +-----------------------------------------------------------------------------*/ + if (randStP->stateP) + free(randStP->stateP); +} + + + |