Use Mersenee Twister instead of libc rand for random numbers

git-svn-id: http://svn.code.sf.net/p/netpbm/code/trunk@4033 9d0c8265-081b-0410-96cb-a4ca84ce46f8

6 files changed, 633 insertions, 0 deletions

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..bea5ea17
--- /dev/null
+++ b/lib/util/randsysrand.c
@@ -0,0 +1,35 @@
+#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..97a1376e
--- /dev/null
+++ b/lib/util/randsysrandom.c
@@ -0,0 +1,34 @@
+#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
+};
+
+