diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/Makefile | 4 | ||||
| -rw-r--r-- | lib/util/Makefile | 4 | ||||
| -rw-r--r-- | lib/util/rand.c | 261 | ||||
| -rw-r--r-- | lib/util/rand.h | 107 | ||||
| -rw-r--r-- | lib/util/randmersenne.c | 192 | ||||
| -rw-r--r-- | lib/util/randsysrand.c | 39 | ||||
| -rw-r--r-- | lib/util/randsysrandom.c | 40 |
7 files changed, 647 insertions, 0 deletions
diff --git a/lib/Makefile b/lib/Makefile index bc758df4..d42658a2 100644 --- a/lib/Makefile +++ b/lib/Makefile @@ -48,6 +48,10 @@ LIBOBJECTS_X = \ util/matrix.o \ util/nsleep.o \ util/nstring.o \ + util/rand.o \ + util/randsysrand.o \ + util/randsysrandom.o \ + util/randmersenne.o \ util/runlength.o \ util/shhopt.o \ util/token.o \ diff --git a/lib/util/Makefile b/lib/util/Makefile index 02119edf..d8e2d135 100644 --- a/lib/util/Makefile +++ b/lib/util/Makefile @@ -19,6 +19,10 @@ UTILOBJECTS = \ matrix.o \ nsleep.o \ nstring.o \ + rand.o \ + randsysrand.o \ + randsysrandom.o \ + randmersenne.o \ runlength.o \ shhopt.o \ token.o \ diff --git a/lib/util/rand.c b/lib/util/rand.c new file mode 100644 index 00000000..2f60de83 --- /dev/null +++ b/lib/util/rand.c @@ -0,0 +1,261 @@ +/* + +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 tests 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; +} + + + +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); +} + + + + diff --git a/lib/util/rand.h b/lib/util/rand.h new file mode 100644 index 00000000..c441890a --- /dev/null +++ b/lib/util/rand.h @@ -0,0 +1,107 @@ +/* Interface header file for random number generator functions in libnetpbm */ + +#ifndef RAND_H_INCLUDED +#define RAND_H_INCLUDED + +#include "netpbm/pm_c_util.h" +#include "netpbm/mallocvar.h" + +/* + Definitions for selecting the random number generator + + The default random number generator is Mersenne Twister. Here we + provide a means to revert to the system rand() generator if need be. + + Glibc provides generators: rand(), random() and drand48(). Each has + its own associated functions. In the Glibc documentation rand() is + called "ISO", random() is called "BSD" and drand48() is called "SVID". + If your system has the glibc documentation installed "info rand" should + get you to the relevant page. The documentation is available online + from: + + https://www.gnu.org/software/libc/manual/html_node/Pseudo_002dRandom-Numbers.html + Pseudo-Random Numbers (The GNU C Library) + + Glibc's choice of name is confusing for what it calls "ISO" rand() + was available in early BSD systems. + + Functions by these names appear on most Unix systems, but generation + formulas and default initial states are known to differ. On systems + which do not use glibc, what is called rand() may have no relation + with the formula of the ISO C standard. Likewise random() may have + no relation with the BSD formula. +*/ + +enum PmRandEngine {PM_RAND_SYS_RAND, /* rand() */ + PM_RAND_SYS_RANDOM, /* random() */ + PM_RAND_SYS_DRAND48, /* drand48() reserved */ + PM_RAND_MERSENNETWISTER /* default */}; + +#ifndef PM_RANDOM_NUMBER_GENERATOR + #define PM_RANDOM_NUMBER_GENERATOR PM_RAND_MERSENNETWISTER +#endif + + +/* Structure to hold random number generator profile and internal state */ + +struct pm_randSt; + +struct pm_rand_vtable { + void + (*init)(struct pm_randSt * const randStP); + + void + (*srand)(struct pm_randSt * const randStP, + unsigned int const seed); + + unsigned long int + (*rand)(struct pm_randSt * const randStP); +}; + +extern struct pm_rand_vtable const pm_randsysrand_vtable; +extern struct pm_rand_vtable const pm_randsysrandom_vtable; +extern struct pm_rand_vtable const pm_randmersenne_vtable; + +struct pm_randSt { + struct pm_rand_vtable vtable; + void * stateP; /* Internal state */ + unsigned int max; + unsigned int seed; + bool gaussCacheValid; + double gaussCache; +}; + +/* Function declarations */ + +extern void +pm_randinit(struct pm_randSt * const randStP); + +extern void +pm_randterm(struct pm_randSt * const randStP); + +extern void +pm_srand(struct pm_randSt * const randStP, + unsigned int const seed); + + +extern void +pm_srand2(struct pm_randSt * const randStP, + bool const withSeed, + unsigned int const seedVal); + +extern unsigned long int +pm_rand(struct pm_randSt * const randStP); + +extern double +pm_drand(struct pm_randSt * const randStP); + +extern void +pm_gaussrand2(struct pm_randSt * const randStP, + double * const r1P, + double * const r2P); + +extern double +pm_gaussrand(struct pm_randSt * const randStP); + + +#endif diff --git a/lib/util/randmersenne.c b/lib/util/randmersenne.c new file mode 100644 index 00000000..34355a23 --- /dev/null +++ b/lib/util/randmersenne.c @@ -0,0 +1,192 @@ +#include "netpbm/pm.h" +#include "netpbm/rand.h" + +/* +++++ Start of Mersenne Twister pseudorandom number generator code +++++ */ + +/* + Original source code from: + http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/VERSIONS/C-LANG/c-lang.html + + A C-program for MT19937, with initialization improved 2002/1/26. + Coded by Takuji Nishimura and Makoto Matsumoto. + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) + + Above conditions apply in the following code to the line which says: + +++++ End of Mersenne Twister pseudorandom number generator code +++++ +*/ + +/* Period parameters */ + +#define MT_N 624 +#define MT_M 397 +#define MT_MATRIX_A 0x9908b0dfUL /* constant vector a */ + +struct MtState { + uint32_t mt[MT_N]; /* the array for the state vector */ + unsigned int mtIndex; +}; + + + +static void +randMtAlloc(struct MtState ** const statePP) { + + struct MtState * stateP; + + MALLOCVAR_NOFAIL(stateP); + + *statePP = stateP; +} + + + +/* 32 bit masks */ + +static uint32_t const FMASK = 0xffffffffUL; /* all bits */ +static uint32_t const UMASK = 0x80000000UL; /* most significant bit */ +static uint32_t const LMASK = 0x7fffffffUL; /* least significant 31 bits */ + + + +static void +srandMt(struct MtState * const stateP, + unsigned int const seed) { +/*----------------------------------------------------------------------------- + Initialize state array mt[MT_N] with seed +-----------------------------------------------------------------------------*/ + unsigned int mtIndex; + uint32_t * const mt = stateP->mt; + + mt[0]= seed & FMASK; + + for (mtIndex = 1; mtIndex < MT_N; ++mtIndex) { + mt[mtIndex] = (1812433253UL * (mt[mtIndex-1] + ^ (mt[mtIndex-1] >> 30)) + mtIndex); + + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + } + + stateP->mtIndex = mtIndex; +} + + + +static unsigned long int +randMt32(struct MtState * const stateP) { +/*---------------------------------------------------------------------------- + Generate a 32 bit random number interval: [0, 0xffffffff] + ----------------------------------------------------------------------------*/ + unsigned int mtIndex; + uint32_t retval; + + if (stateP->mtIndex >= MT_N) { + /* generate N words at one time */ + uint32_t * const mt = stateP->mt; + uint32_t const mag01[2]={0x0UL, MT_MATRIX_A}; + /* mag01[x] = x * MT_MATRIX_A for x=0, 1 */ + + int k; + uint32_t y; + + if (stateP->mtIndex >= MT_N+1) { + pm_error("Internal error in Mersenne Twister random number" + "generator"); + } + + for (k = 0; k < MT_N-MT_M; ++k) { + y = (mt[k] & UMASK) | (mt[k+1] & LMASK); + mt[k] = mt[k + MT_M] ^ (y >> 1) ^ mag01[y & 0x1UL]; + } + for (; k < MT_N-1; ++k) { + y = (mt[k] & UMASK) | (mt[k+1] & LMASK); + mt[k] = mt[k+(MT_M-MT_N)] ^ (y >> 1) ^ mag01[y & 0x1UL]; + } + y = (mt[MT_N - 1] & UMASK) | (mt[0] & LMASK); + mt[MT_N - 1] = mt[MT_M - 1] ^ (y >> 1) ^ mag01[y & 0x1UL]; + + mtIndex = 0; + } else + mtIndex = stateP->mtIndex; + + retval = stateP->mt[mtIndex]; + + /* Tempering */ + retval ^= (retval >> 11); + retval ^= (retval << 7) & 0x9d2c5680UL; + retval ^= (retval << 15) & 0xefc60000UL; + retval ^= (retval >> 18); + + stateP->mtIndex = mtIndex + 1; + + return retval; +} + +/* +++++ End of Mersenne Twister pseudorandom number generator code +++++ */ + + +static void +vinit(struct pm_randSt * const randStP) { + + randMtAlloc((struct MtState ** const) &randStP->stateP); + randStP->max = 0xffffffffUL; +} + + + +static void +vsrand(struct pm_randSt * const randStP, + unsigned int const seed) { + + srandMt(randStP->stateP, seed); +} + + + +static unsigned long int +vrand(struct pm_randSt * const randStP) { + + return randMt32(randStP->stateP); +} + + + +struct pm_rand_vtable const pm_randmersenne_vtable = { + &vinit, + &vsrand, + &vrand +}; + + diff --git a/lib/util/randsysrand.c b/lib/util/randsysrand.c new file mode 100644 index 00000000..f97a5d3c --- /dev/null +++ b/lib/util/randsysrand.c @@ -0,0 +1,39 @@ +#include <stdlib.h> + +#include "netpbm/rand.h" + + + +static void +vinit(struct pm_randSt * const randStP) { + + randStP->max = RAND_MAX; + randStP->stateP = NULL; +} + + + +static void +vsrand(struct pm_randSt * const randStP, + unsigned int const seed) { + + srand(seed); +} + + + +static unsigned long int +vrand(struct pm_randSt * const randStP) { + + return rand(); +} + + + +struct pm_rand_vtable const pm_randsysrand_vtable = { + &vinit, + &vsrand, + &vrand +}; + + diff --git a/lib/util/randsysrandom.c b/lib/util/randsysrandom.c new file mode 100644 index 00000000..3e38b0ae --- /dev/null +++ b/lib/util/randsysrandom.c @@ -0,0 +1,40 @@ +#define _DEFAULT_SOURCE /* New name for SVID & BSD source defines */ +#define _XOPEN_SOURCE 500 /* Make sure random() is in stdlib.h */ +#define _BSD_SOURCE /* Make sure random() is in stdlib.h */ + +#include <stdlib.h> + +#include "netpbm/rand.h" + +static void +vinit(struct pm_randSt * const randStP) { + + randStP->max = RAND_MAX; + randStP->stateP = NULL; +} + + + +static void +vsrand(struct pm_randSt * const randStP, + unsigned int const seed) { + + srandom(seed); +} + + + +static unsigned long int +vrand(struct pm_randSt * const randStP) { + + return random(); +} + + +struct pm_rand_vtable const pm_randsysrandom_vtable = { + &vinit, + &vsrand, + &vrand +}; + + |