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

25 files changed, 1828 insertions, 797 deletions

diff --git a/common.mk b/common.mk
index 749488c2..d2fa50d8 100644
--- a/common.mk
+++ b/common.mk
@@ -150,7 +150,7 @@ IMPORTINC_LIB_HEADERS := \
 IMPORTINC_LIB_UTIL_HEADERS := \
   bitarith.h bitio.h bitreverse.h filename.h intcode.h floatcode.h io.h \
   matrix.h mallocvar.h \
-  nsleep.h nstring.h pm_c_util.h runlength.h shhopt.h token.h
+  nsleep.h nstring.h pm_c_util.h rand.h runlength.h shhopt.h token.h
 
 IMPORTINC_HEADERS := \
   $(IMPORTINC_ROOT_HEADERS) \
diff --git a/converter/other/pgmtopbm.c b/converter/other/pgmtopbm.c
index 64dc814b..d5f67a06 100644
--- a/converter/other/pgmtopbm.c
+++ b/converter/other/pgmtopbm.c
@@ -17,6 +17,7 @@
 #include "pgm.h"
 #include "dithers.h"
 #include "mallocvar.h"
+#include "rand.h"
 
 enum halftone {QT_FS, QT_THRESH, QT_DITHER8, QT_CLUSTER, QT_HILBERT};
 
@@ -462,14 +463,19 @@ createFsConverter(unsigned int const cols,
     /* Initialize Floyd-Steinberg error vectors. */
     MALLOCARRAY_NOFAIL(stateP->thiserr, cols + 2);
     MALLOCARRAY_NOFAIL(stateP->nexterr, cols + 2);
-    srand(randomSeedSpec ? randomSeed : pm_randseed());
 
     {
         /* (random errors in [-fs_scale/8 .. fs_scale/8]) */
         unsigned int col;
+        struct pm_randSt randSt;
+        pm_randinit(&randSt);
+        pm_srand2(&randSt, randomSeedSpec, randomSeed);
+
         for (col = 0; col < cols + 2; ++col)
             stateP->thiserr[col] =
-                (long)(rand() % fs_scale - half_fs_scale) / 4;
+                (long)(pm_rand(&randSt) % fs_scale - half_fs_scale) / 4;
+
+        pm_randterm(&randSt);
     }
 
     stateP->fs_forward = TRUE;
diff --git a/editor/pamaddnoise.c b/editor/pamaddnoise.c
index 20c68d99..9ca80394 100644
--- a/editor/pamaddnoise.c
+++ b/editor/pamaddnoise.c
@@ -33,18 +33,20 @@
 
 #include "pm_c_util.h"
 #include "mallocvar.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "pm_gamma.h"
 #include "pam.h"
 
 static double const EPSILON = 1.0e-5;
+static double const SALT_RATIO = 0.5;
 
 
 
 static double
-rand1() {
+rand1(struct pm_randSt * const randStP) {
 
-    return (double)rand()/RAND_MAX;
+    return (double)pm_rand(randStP)/RAND_MAX;
 }
 
 
@@ -67,6 +69,7 @@ struct CmdlineInfo {
 
     enum NoiseType noiseType;
 
+    unsigned int seedSpec;
     unsigned int seed;
 
     float lambda;
@@ -119,7 +122,7 @@ parseCommandLine(int argc, const char ** const argv,
 
     unsigned int option_def_index;
 
-    unsigned int typeSpec, seedSpec, lambdaSpec, lsigmaSpec, mgsigmaSpec,
+    unsigned int typeSpec, lambdaSpec, lsigmaSpec, mgsigmaSpec,
         sigma1Spec, sigma2Spec, toleranceSpec;
 
     const char * type;
@@ -128,21 +131,21 @@ parseCommandLine(int argc, const char ** const argv,
 
     option_def_index = 0;   /* incremented by OPTENT3 */
     OPTENT3(0,   "type",            OPT_STRING,   &type,
-            &typeSpec,         0);
+            &typeSpec,           0);
     OPTENT3(0,   "seed",            OPT_UINT,     &cmdlineP->seed,
-            &seedSpec,         0);
+            &cmdlineP->seedSpec, 0);
     OPTENT3(0,   "lambda",          OPT_FLOAT,    &cmdlineP->lambda,
-            &lambdaSpec,       0);
+            &lambdaSpec,         0);
     OPTENT3(0,   "lsigma",          OPT_FLOAT,    &cmdlineP->lsigma,
-            &lsigmaSpec,       0);
+            &lsigmaSpec,         0);
     OPTENT3(0,   "mgsigma",         OPT_FLOAT,    &cmdlineP->mgsigma,
-            &mgsigmaSpec,      0);
+            &mgsigmaSpec,        0);
     OPTENT3(0,   "sigma1",          OPT_FLOAT,    &cmdlineP->sigma1,
-            &sigma1Spec,       0);
+            &sigma1Spec,         0);
     OPTENT3(0,   "sigma2",          OPT_FLOAT,    &cmdlineP->sigma2,
-            &sigma2Spec,       0);
+            &sigma2Spec,         0);
     OPTENT3(0,   "tolerance",       OPT_FLOAT,    &cmdlineP->tolerance,
-            &toleranceSpec,    0);
+            &toleranceSpec,      0);
 
     opt.opt_table = option_def;
     opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
@@ -193,7 +196,7 @@ parseCommandLine(int argc, const char ** const argv,
     if (!toleranceSpec)
         cmdlineP->tolerance = 0.10;
 
-    if (!seedSpec)
+    if (!cmdlineP->seedSpec)
         cmdlineP->seed = pm_randseed();
 
     if (argc-1 > 1)
@@ -211,11 +214,12 @@ parseCommandLine(int argc, const char ** const argv,
 
 
 static void
-addGaussianNoise(sample   const maxval,
-                 sample   const origSample,
-                 sample * const newSampleP,
-                 float    const sigma1,
-                 float    const sigma2) {
+addGaussianNoise(sample             const maxval,
+                 sample             const origSample,
+                 sample *           const newSampleP,
+                 float              const sigma1,
+                 float              const sigma2,
+                 struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    Add Gaussian noise.
 
@@ -225,11 +229,11 @@ addGaussianNoise(sample   const maxval,
     double x1, x2, xn, yn;
     double rawNewSample;
 
-    x1 = rand1();
+    x1 = rand1(randStP);
 
     if (x1 == 0.0)
         x1 = 1.0;
-    x2 = rand1();
+    x2 = rand1(randStP);
     xn = sqrt(-2.0 * log(x1)) * cos(2.0 * M_PI * x2);
     yn = sqrt(-2.0 * log(x1)) * sin(2.0 * M_PI * x2);
 
@@ -242,40 +246,42 @@ addGaussianNoise(sample   const maxval,
 
 
 static void
-addImpulseNoise(sample   const maxval,
-                sample   const origSample,
-                sample * const newSampleP,
-                float    const tolerance) {
+addImpulseNoise(sample             const maxval,
+                sample             const origSample,
+                sample *           const newSampleP,
+                float              const tolerance,
+                double             const saltRatio,
+                struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    Add impulse (salt and pepper) noise
 -----------------------------------------------------------------------------*/
 
-    double const low_tol  = tolerance / 2.0;
-    double const high_tol = 1.0 - (tolerance / 2.0);
-    double const sap = rand1();
+    double const pepperRatio = 1.0 - saltRatio;
+    double const loTolerance = tolerance * pepperRatio;
+    double const hiTolerance = 1.0 - tolerance * saltRatio;
+    double const sap         = rand1(randStP);
 
-    if (sap < low_tol)
-        *newSampleP = 0;
-    else if ( sap >= high_tol )
-        *newSampleP = maxval;
-    else
-        *newSampleP = origSample;
+    *newSampleP =
+        sap < loTolerance ? 0 :
+        sap >= hiTolerance? maxval :
+        origSample;
 }
 
 
 
 static void
-addLaplacianNoise(sample   const maxval,
-                  double   const infinity,
-                  sample   const origSample,
-                  sample * const newSampleP,
-                  float    const lsigma) {
+addLaplacianNoise(sample             const maxval,
+                  double             const infinity,
+                  sample             const origSample,
+                  sample *           const newSampleP,
+                  float              const lsigma,
+                  struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    Add Laplacian noise
 
    From Pitas' book.
 -----------------------------------------------------------------------------*/
-    double const u = rand1();
+    double const u = rand1(randStP);
 
     double rawNewSample;
 
@@ -297,11 +303,12 @@ addLaplacianNoise(sample   const maxval,
 
 
 static void
-addMultiplicativeGaussianNoise(sample   const maxval,
-                               double   const infinity,
-                               sample   const origSample,
-                               sample * const newSampleP,
-                               float    const mgsigma) {
+addMultiplicativeGaussianNoise(sample             const maxval,
+                               double             const infinity,
+                               sample             const origSample,
+                               sample *           const newSampleP,
+                               float              const mgsigma,
+                               struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    Add multiplicative Gaussian noise
 
@@ -311,14 +318,14 @@ addMultiplicativeGaussianNoise(sample   const maxval,
     double rawNewSample;
 
     {
-        double const uniform = rand1();
+        double const uniform = rand1(randStP);
         if (uniform <= EPSILON)
             rayleigh = infinity;
         else
             rayleigh = sqrt(-2.0 * log( uniform));
     }
     {
-        double const uniform = rand1();
+        double const uniform = rand1(randStP);
         gauss = rayleigh * cos(2.0 * M_PI * uniform);
     }
     rawNewSample = origSample + (origSample * mgsigma * gauss);
@@ -363,10 +370,11 @@ poissonPmf(double       const lambda,
 
 
 static void
-addPoissonNoise(struct pam * const pamP,
-                sample       const origSample,
-                sample *     const newSampleP,
-                float        const lambdaOfMaxval) {
+addPoissonNoise(struct pam *       const pamP,
+                sample             const origSample,
+                sample *           const newSampleP,
+                float              const lambdaOfMaxval,
+                struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    Add Poisson noise
 -----------------------------------------------------------------------------*/
@@ -376,7 +384,7 @@ addPoissonNoise(struct pam * const pamP,
 
     double const lambda  = origSampleIntensity * lambdaOfMaxval;
 
-    double const u = rand1();
+    double const u = rand1(randStP);
 
     /* We now apply the inverse CDF (cumulative distribution function) of the
        Poisson distribution to uniform random variable 'u' to get a Poisson
@@ -416,12 +424,14 @@ main(int argc, const char ** argv) {
     const tuple * newtuplerow;
     unsigned int row;
     double infinity;
+    struct pm_randSt randSt;
 
     pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
 
-    srand(cmdline.seed);
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, cmdline.seedSpec, cmdline.seed);
 
     ifP = pm_openr(cmdline.inputFileName);
 
@@ -448,35 +458,40 @@ main(int argc, const char ** argv) {
                     addGaussianNoise(inpam.maxval,
                                      tuplerow[col][plane],
                                      &newtuplerow[col][plane],
-                                     cmdline.sigma1, cmdline.sigma2);
+                                     cmdline.sigma1, cmdline.sigma2,
+                                     &randSt);
                     break;
 
                 case NOISETYPE_IMPULSE:
                     addImpulseNoise(inpam.maxval,
                                     tuplerow[col][plane],
                                     &newtuplerow[col][plane],
-                                    cmdline.tolerance);
+                                    cmdline.tolerance, SALT_RATIO,
+                                    &randSt);
                    break;
 
                 case NOISETYPE_LAPLACIAN:
                     addLaplacianNoise(inpam.maxval, infinity,
                                       tuplerow[col][plane],
                                       &newtuplerow[col][plane],
-                                      cmdline.lsigma);
+                                      cmdline.lsigma,
+                                      &randSt);
                     break;
 
                 case NOISETYPE_MULTIPLICATIVE_GAUSSIAN:
                     addMultiplicativeGaussianNoise(inpam.maxval, infinity,
                                                    tuplerow[col][plane],
                                                    &newtuplerow[col][plane],
-                                                   cmdline.mgsigma);
+                                                   cmdline.mgsigma,
+                                                   &randSt);
                     break;
 
                 case NOISETYPE_POISSON:
                     addPoissonNoise(&inpam,
                                     tuplerow[col][plane],
                                     &newtuplerow[col][plane],
-                                    cmdline.lambda);
+                                    cmdline.lambda,
+                                    &randSt);
                     break;
 
                 }
@@ -484,6 +499,7 @@ main(int argc, const char ** argv) {
         }
         pnm_writepamrow(&outpam, newtuplerow);
     }
+    pm_randterm(&randSt);
     pnm_freepamrow(newtuplerow);
     pnm_freepamrow(tuplerow);
 
diff --git a/editor/pamditherbw.c b/editor/pamditherbw.c
index ae91a26f..694b2c21 100644
--- a/editor/pamditherbw.c
+++ b/editor/pamditherbw.c
@@ -14,12 +14,14 @@
 #include <string.h>
 
 #include "pm_c_util.h"
-#include "pam.h"
-#include "dithers.h"
+#include "rand.h"
 #include "mallocvar.h"
 #include "shhopt.h"
+#include "pam.h"
+#include "dithers.h"
 #include "pm_gamma.h"
 
+
 enum halftone {QT_FS,
                QT_ATKINSON,
                QT_THRESH,
@@ -588,7 +590,9 @@ fsDestroy(struct converter * const converterP) {
 
 static struct converter
 createFsConverter(struct pam * const graypamP,
-                  float        const threshFraction) {
+                  float        const threshFraction,
+                  bool         const randomseedSpec,
+                  unsigned int const randomseed) {
 
     struct fsState * stateP;
     struct converter converter;
@@ -605,9 +609,17 @@ createFsConverter(struct pam * const graypamP,
 
     {
         /* (random errors in [-1/8 .. 1/8]) */
+
         unsigned int col;
+        struct pm_randSt randSt;
+
+        pm_randinit(&randSt);
+        pm_srand2(&randSt, randomseedSpec, randomseed);
+
         for (col = 0; col < graypamP->width + 2; ++col)
-            stateP->thiserr[col] = ((float)rand()/RAND_MAX - 0.5) / 4;
+            stateP->thiserr[col] = (pm_drand(&randSt) - 0.5) / 4;
+
+        pm_randterm(&randSt);
     }
 
     stateP->halfWhite = threshFraction;
@@ -725,7 +737,9 @@ atkinsonDestroy(struct converter * const converterP) {
 
 static struct converter
 createAtkinsonConverter(struct pam * const graypamP,
-                        float        const threshFraction) {
+                        float        const threshFraction,
+                        bool         const randomseedSpec,
+                        unsigned int const randomseed) {
 
     struct atkinsonState * stateP;
     struct converter converter;
@@ -743,11 +757,18 @@ createAtkinsonConverter(struct pam * const graypamP,
     {
         /* (random errors in [-1/8 .. 1/8]) */
         unsigned int col;
+        struct pm_randSt randSt;
+
+        pm_randinit(&randSt);
+        pm_srand2(&randSt, randomseedSpec, randomseed);
+
         for (col = 0; col < graypamP->width + 2; ++col) {
-            stateP->error[0][col] = ((float)rand()/RAND_MAX - 0.5) / 4;
+            stateP->error[0][col] = (pm_drand(&randSt) - 0.5) / 4;
             stateP->error[1][col] = 0.0;
             stateP->error[2][col] = 0.0;
         }
+
+        pm_randterm(&randSt);
     }
 
     stateP->halfWhite = threshFraction;
@@ -934,8 +955,6 @@ main(int argc, char *argv[]) {
 
     parseCommandLine(argc, argv, &cmdline);
 
-    srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
-
     ifP = pm_openr(cmdline.inputFilespec);
 
     if (cmdline.halftone == QT_HILBERT)
@@ -956,10 +975,14 @@ main(int argc, char *argv[]) {
 
         switch (cmdline.halftone) {
         case QT_FS:
-            converter = createFsConverter(&graypam, cmdline.threshval);
+            converter = createFsConverter(&graypam, cmdline.threshval,
+                                          cmdline.randomseedSpec,
+                                          cmdline.randomseed);
             break;
         case QT_ATKINSON:
-            converter = createAtkinsonConverter(&graypam, cmdline.threshval);
+            converter = createAtkinsonConverter(&graypam, cmdline.threshval,
+                                                cmdline.randomseedSpec,
+                                                cmdline.randomseed);
             break;
         case QT_THRESH:
             converter = createThreshConverter(&graypam, cmdline.threshval);
diff --git a/editor/pammixmulti.c b/editor/pammixmulti.c
index f5012d7a..2b45d807 100644
--- a/editor/pammixmulti.c
+++ b/editor/pammixmulti.c
@@ -13,6 +13,7 @@
 #include "shhopt.h"
 #include "mallocvar.h"
 #include "nstring.h"
+#include "rand.h"
 
 typedef enum {
     BLEND_AVERAGE,   /* Take the average color of all pixels */
@@ -36,6 +37,8 @@ struct ProgramState {
         /* Standard deviation when selecting images via a mask */
     unsigned long ** imageWeights;
         /* Per-image weights as a function of grayscale level */
+    struct pm_randSt randSt;
+        /* Random number generator parameters and internal state */
 };
 
 
@@ -134,9 +137,9 @@ parseCommandLine(int argc, const char ** argv,
 }
 
 static void
-openInputFiles(unsigned int          const inFileCt,
-               const char **         const inFileName,
-               struct ProgramState * const stateP) {
+initInput(unsigned int          const inFileCt,
+          const char **         const inFileName,
+          struct ProgramState * const stateP) {
 /*----------------------------------------------------------------------------
   Open all of the input files.
 
@@ -180,6 +183,25 @@ openInputFiles(unsigned int          const inFileCt,
 }
 
 
+
+static void
+termInput(struct ProgramState * const stateP) {
+/*----------------------------------------------------------------------------
+  Deallocate all of the resources we allocated.
+-----------------------------------------------------------------------------*/
+    unsigned int i;
+
+    for (i = 0; i < stateP->inFileCt; ++i) {
+        pnm_freepamrow(stateP->inTupleRows[i]);
+        pm_close(stateP->inPam[i].file);
+    }
+
+    free(stateP->inTupleRows);
+    free(stateP->inPam);
+}
+
+
+
 static void
 initMask(const char *          const maskFileName,
          struct ProgramState * const stateP) {
@@ -233,6 +255,15 @@ initOutput(FILE *                const ofP,
 }
 
 
+static void
+termOutput(struct ProgramState * const stateP) {
+
+    free(stateP->outTupleRow);
+
+    pm_close(stateP->outPam.file);
+}
+
+
 
 static void
 blendTuplesRandom(struct ProgramState * const stateP,
@@ -243,8 +274,8 @@ blendTuplesRandom(struct ProgramState * const stateP,
   from a random input image.
 -----------------------------------------------------------------------------*/
     unsigned int const depth = stateP->inPam[0].depth;
-    unsigned int const img = (unsigned int) (rand() % stateP->inFileCt);
-
+    unsigned int const img = (unsigned int) (pm_rand(&stateP->randSt) %
+                                             stateP->inFileCt);
     unsigned int samp;
 
     for (samp = 0; samp < depth; ++samp)
@@ -276,23 +307,26 @@ blendTuplesAverage(struct ProgramState * const stateP,
 
 
 
+#if 0
 static void
-randomNormal2(double * const r1P,
-              double * const r2P) {
+randomNormal2(double *           const r1P,
+              double *           const r2P,
+              struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
   Return two normally distributed random numbers.
 -----------------------------------------------------------------------------*/
     double u1, u2;
 
     do {
-        u1 = drand48();
-        u2 = drand48();
+        u1 = drand48(randStP);
+        u2 = drand48(randStP);
     }
     while (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);
 }
+#endif
 
 
 
@@ -332,7 +366,8 @@ precomputeImageWeights(struct ProgramState * const stateP,
             double r[2];
             unsigned int k;
 
-            randomNormal2(&r[0], &r[1]);
+            pm_gaussrand2(&stateP->randSt, &r[0], &r[1]);
+
             for (k = 0; k < 2; ++k) {
                 int const img =
                     r[k] * sigma + pctGray * stateP->inFileCt * 0.999999;
@@ -453,26 +488,6 @@ blendImages(BlendType             const blend,
 
 
 
-static void
-termState(struct ProgramState * const stateP) {
-/*----------------------------------------------------------------------------
-  Deallocate all of the resources we allocated.
------------------------------------------------------------------------------*/
-    unsigned int i;
-
-    for (i = 0; i < stateP->inFileCt; ++i) {
-        pnm_freepamrow(stateP->inTupleRows[i]);
-        pm_close(stateP->inPam[i].file);
-    }
-
-    free(stateP->outTupleRow);
-    free(stateP->inTupleRows);
-    free(stateP->inPam);
-    pm_close(stateP->outPam.file);
-}
-
-
-
 int
 main(int argc, const char * argv[]) {
 
@@ -483,13 +498,14 @@ main(int argc, const char * argv[]) {
 
     parseCommandLine(argc, argv, &cmdline);
 
-    srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
-
-    openInputFiles(cmdline.inFileNameCt, cmdline.inFileName, &state);
+    initInput(cmdline.inFileNameCt, cmdline.inFileName, &state);
 
     if (cmdline.blend == BLEND_MASK)
         initMask(cmdline.maskfile, &state);
 
+    pm_randinit(&state.randSt);
+    pm_srand2(&state.randSt, cmdline.randomseedSpec, cmdline.randomseed);
+
     initOutput(stdout, &state);
 
     if (cmdline.blend == BLEND_MASK)
@@ -497,10 +513,14 @@ main(int argc, const char * argv[]) {
 
     blendImages(cmdline.blend, &state);
 
+    termOutput(&state);
+
+    pm_randterm(&state.randSt);
+
     if (cmdline.blend == BLEND_MASK)
         termMask(&state);
 
-    termState(&state);
+    termInput(&state);
 
     freeCmdline(&cmdline);
 
diff --git a/editor/pamrecolor.c b/editor/pamrecolor.c
index 8c5bce12..644d022d 100644
--- a/editor/pamrecolor.c
+++ b/editor/pamrecolor.c
@@ -32,6 +32,7 @@
 
 #include "mallocvar.h"
 #include "nstring.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "pam.h"
 
@@ -175,17 +176,25 @@ explicitlyColorRow(struct pam *   const pamP,
 
 static void
 randomlyColorRow(struct pam *   const pamP,
-                 tuplen *       const rowData) {
+                 tuplen *       const rowData,
+                 bool           const randomseedSpec,
+                 unsigned int   const randomseed) {
 /*----------------------------------------------------------------------
   Assign each tuple in a row a random color.
 ------------------------------------------------------------------------*/
     unsigned int col;
+    struct pm_randSt randSt;
+
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, randomseedSpec, randomseed);
 
     for (col = 0; col < pamP->width; ++col) {
-        rowData[col][PAM_RED_PLANE] = rand() / (float)RAND_MAX;
-        rowData[col][PAM_GRN_PLANE] = rand() / (float)RAND_MAX;
-        rowData[col][PAM_BLU_PLANE] = rand() / (float)RAND_MAX;
+        rowData[col][PAM_RED_PLANE] = pm_drand(&randSt);
+        rowData[col][PAM_GRN_PLANE] = pm_drand(&randSt);
+        rowData[col][PAM_BLU_PLANE] = pm_drand(&randSt);
     }
+
+    pm_randterm(&randSt);
 }
 
 
@@ -462,8 +471,6 @@ main(int argc, const char *argv[]) {
 
     parseCommandLine(argc, argv, &cmdline);
 
-    srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
-
     ifP = pm_openr(cmdline.inputFileName);
     inPam.comment_p = &comments;
     pnm_readpaminit(ifP, &inPam, PAM_STRUCT_SIZE(comment_p));
@@ -505,7 +512,8 @@ main(int argc, const char *argv[]) {
             if (cmdline.targetcolorSpec)
                 explicitlyColorRow(&colorPam, colorRow, cmdline.targetcolor);
             else
-                randomlyColorRow(&colorPam, colorRow);
+                randomlyColorRow(&colorPam, colorRow,
+                                 cmdline.randomseedSpec, cmdline.randomseed);
         }
         recolorRow(&inPam, inRow,
                    &cmdline.color2gray, colorRow,
diff --git a/editor/pamrubber.c b/editor/pamrubber.c
index 602701ec..fda31203 100644
--- a/editor/pamrubber.c
+++ b/editor/pamrubber.c
@@ -1,20 +1,18 @@
-/*----------------------------------------------------------------------------*/
-
-/* pamrubber.c - transform images using Rubber Sheeting algorithm
-**               see: http://www.schaik.com/netpbm/rubber/
-**
-** Copyright (C) 2011 by Willem van Schaik (willem@schaik.com)
-**
-** Permission to use, copy, modify, and distribute this software and its
-** documentation for any purpose and without fee is hereby granted, provided
-** that the above copyright notice appear in all copies and that both that
-** copyright notice and this permission notice appear in supporting
-** documentation.  This software is provided "as is" without express or
-** implied warranty.
-*/
-
-/*----------------------------------------------------------------------------*/
-
+/*=============================================================================
+                              pamrubber
+===============================================================================
+  Transform images using Rubber Sheeting algorithm
+  See: http://www.schaik.com/netpbm/rubber/
+
+  Copyright (C) 2011 by Willem van Schaik (willem@schaik.com)
+
+  Permission to use, copy, modify, and distribute this software and its
+  documentation for any purpose and without fee is hereby granted, provided
+  that the above copyright notice appear in all copies and that both that
+  copyright notice and this permission notice appear in supporting
+  documentation.  This software is provided "as is" without express or
+  implied warranty.
+=============================================================================*/
 #include <stdio.h>
 #include <stdlib.h>
 #include <unistd.h>
@@ -25,12 +23,12 @@
 
 #include "pm_c_util.h"
 #include "mallocvar.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "pam.h"
 #include "pamdraw.h"
 
 
-
 typedef struct {
   double x;
   double y;
@@ -54,7 +52,7 @@ typedef struct {
     point br;  /* bottom right */
 } quadrilateral;
 
-struct cmdlineInfo {
+struct CmdlineInfo {
     unsigned int nCP;
     point        oldCP[4];
     point        newCP[4];
@@ -64,14 +62,14 @@ struct cmdlineInfo {
     unsigned int frame;
     unsigned int linear;
     unsigned int verbose;
-    unsigned int randseedSpec;
-    unsigned int randseed;
+    unsigned int randomseedSpec;
+    unsigned int randomseed;
 };
 
 
 static void
 parseCmdline(int argc, const char ** argv,
-             struct cmdlineInfo * const cmdlineP) {
+             struct CmdlineInfo * const cmdlineP) {
 
 /* parse all parameters from the command line */
 
@@ -92,10 +90,10 @@ parseCmdline(int argc, const char ** argv,
     OPTENT3(0, "frame",    OPT_FLAG, NULL, &cmdlineP->frame,    0);
     OPTENT3(0, "linear",   OPT_FLAG, NULL, &cmdlineP->linear,   0);
     OPTENT3(0, "verbose",  OPT_FLAG, NULL, &cmdlineP->verbose,  0);
-    OPTENT3(0, "randseed", OPT_UINT, &cmdlineP->randseed,
-            &cmdlineP->randseedSpec, 0);
-    OPTENT3(0, "randomseed", OPT_UINT, &cmdlineP->randseed,
-            &cmdlineP->randseedSpec, 0);
+    OPTENT3(0, "randseed", OPT_UINT, &cmdlineP->randomseed,
+            &cmdlineP->randomseedSpec, 0);
+    OPTENT3(0, "randomseed", OPT_UINT, &cmdlineP->randomseed,
+            &cmdlineP->randomseedSpec, 0);
 
     opt.opt_table = option_def;
     opt.short_allowed = FALSE;  /* we have no short (old-fashioned) options */
@@ -339,28 +337,29 @@ windtriangle(triangle * const tP,
 
 
 static double
-tiny(void) {
+tiny(struct pm_randSt * const randStP) {
 
-    if (rand() % 2)
-        return +1E-6 * (double) ((rand() % 90) + 9);
+    if (pm_rand(randStP) % 2)
+        return +1E-6 * (double) ((pm_rand(randStP) % 90) + 9);
     else
-        return -1E-6 * (double) ((rand() % 90) + 9);
+        return -1E-6 * (double) ((pm_rand(randStP) % 90) + 9);
 }
 
 
 
 static void
 angle(point * const p1P,
-      point * const p2P) {
+      point * const p2P,
+      struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    Move *p2P slightly if necessary to make sure the line (*p1P, *p2P)
    is not horizontal or vertical.
 -----------------------------------------------------------------------------*/
     if (p1P->x == p2P->x) { /* vertical line */
-        p2P->x += tiny();
+        p2P->x += tiny(randStP);
     }
     if (p1P->y == p2P->y) { /* horizontal line */
-        p2P->y += tiny();
+        p2P->y += tiny(randStP);
     }
 }
 
@@ -720,8 +719,10 @@ static void drawClippedTriangle(const struct pam * const pamP,
 
 
 static void
-prepTrig(int const wd,
-         int const ht) {
+prepTrig(int          const wd,
+         int          const ht,
+         bool         const randomseedSpec,
+         unsigned int const randomseed) {
 
 /* create triangles using control points */
 
@@ -731,16 +732,20 @@ prepTrig(int const wd,
     point c2p1, c2p2, c2p3, c2p4;
     line l1, l2;
     point p0;
+    struct pm_randSt randSt;
+
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, randomseedSpec, randomseed);
 
-    rtl1 = makepoint(0.0 + tiny(),               0.0 + tiny());
-    rtr1 = makepoint((double) wd - 1.0 + tiny(), 0.0 + tiny());
-    rbl1 = makepoint(0.0 + tiny(),               (double) ht - 1.0 + tiny());
-    rbr1 = makepoint((double) wd - 1.0 + tiny(), (double) ht - 1.0 + tiny());
+    rtl1 = makepoint(0.0 + tiny(&randSt),               0.0 + tiny(&randSt));
+    rtr1 = makepoint((double) wd - 1.0 + tiny(&randSt), 0.0 + tiny(&randSt));
+    rbl1 = makepoint(0.0 + tiny(&randSt),               (double) ht - 1.0 + tiny(&randSt));
+    rbr1 = makepoint((double) wd - 1.0 + tiny(&randSt), (double) ht - 1.0 + tiny(&randSt));
 
-    rtl2 = makepoint(0.0 + tiny(),               0.0 + tiny());
-    rtr2 = makepoint((double) wd - 1.0 + tiny(), 0.0 + tiny());
-    rbl2 = makepoint(0.0 + tiny(),               (double) ht - 1.0 + tiny());
-    rbr2 = makepoint((double) wd - 1.0 + tiny(), (double) ht - 1.0 + tiny());
+    rtl2 = makepoint(0.0 + tiny(&randSt),               0.0 + tiny(&randSt));
+    rtr2 = makepoint((double) wd - 1.0 + tiny(&randSt), 0.0 + tiny(&randSt));
+    rbl2 = makepoint(0.0 + tiny(&randSt),               (double) ht - 1.0 + tiny(&randSt));
+    rbr2 = makepoint((double) wd - 1.0 + tiny(&randSt), (double) ht - 1.0 + tiny(&randSt));
 
     if (nCP == 1) {
         c1p1 = oldCP[0];
@@ -772,8 +777,8 @@ prepTrig(int const wd,
         c2p2 = newCP[1];
 
         /* check for hor/ver edges */
-        angle (&c1p1, &c1p2);
-        angle (&c2p1, &c2p2);
+        angle (&c1p1, &c1p2, &randSt);
+        angle (&c2p1, &c2p2, &randSt);
 
         /* connect two control points to corners to get 6 triangles */
         /* left side */
@@ -811,13 +816,13 @@ prepTrig(int const wd,
         /* Move vertices slightly if necessary to make sure no edge is
            horizontal or vertical.
         */
-        angle(&c1p1, &c1p2);
-        angle(&c1p2, &c1p3);
-        angle(&c1p3, &c1p1);
+        angle(&c1p1, &c1p2, &randSt);
+        angle(&c1p2, &c1p3, &randSt);
+        angle(&c1p3, &c1p1, &randSt);
 
-        angle(&c2p1, &c2p2);
-        angle(&c2p2, &c2p3);
-        angle(&c2p3, &c2p1);
+        angle(&c2p1, &c2p2, &randSt);
+        angle(&c2p2, &c2p3, &randSt);
+        angle(&c2p3, &c2p1, &randSt);
 
         if (windtriangle(&tri1s[0], c1p1, c1p2, c1p3)) {
             tri2s[0] = maketriangle(c2p1, c2p2, c2p3);
@@ -871,19 +876,19 @@ prepTrig(int const wd,
         c2p4 = newCP[3];
 
         /* check for hor/ver edges */
-        angle (&c1p1, &c1p2);
-        angle (&c1p2, &c1p3);
-        angle (&c1p3, &c1p4);
-        angle (&c1p4, &c1p1);
-        angle (&c1p1, &c1p3);
-        angle (&c1p2, &c1p4);
-
-        angle (&c2p1, &c2p2);
-        angle (&c2p2, &c2p3);
-        angle (&c2p3, &c2p4);
-        angle (&c2p4, &c2p1);
-        angle (&c2p1, &c2p3);
-        angle (&c2p2, &c2p4);
+        angle (&c1p1, &c1p2, &randSt);
+        angle (&c1p2, &c1p3, &randSt);
+        angle (&c1p3, &c1p4, &randSt);
+        angle (&c1p4, &c1p1, &randSt);
+        angle (&c1p1, &c1p3, &randSt);
+        angle (&c1p2, &c1p4, &randSt);
+
+        angle (&c2p1, &c2p2, &randSt);
+        angle (&c2p2, &c2p3, &randSt);
+        angle (&c2p3, &c2p4, &randSt);
+        angle (&c2p4, &c2p1, &randSt);
+        angle (&c2p1, &c2p3, &randSt);
+        angle (&c2p2, &c2p4, &randSt);
 
         /*-------------------------------------------------------------------*/
         /*        -1-      -2-        -3-      -4-        -5-      -6-       */
@@ -979,6 +984,8 @@ prepTrig(int const wd,
                      &tri2s[9], c2p3, c2p1, rtl2, rtr2, rbl2, rbr2);
         nTri = 10;
     }
+
+    pm_randterm(&randSt);
 }
 
 
@@ -1277,7 +1284,7 @@ warpQuad(point   const p2,
 
 
 static void
-setGlobalCP(struct cmdlineInfo const cmdline) {
+setGlobalCP(struct CmdlineInfo const cmdline) {
 
     unsigned int i;
 
@@ -1392,7 +1399,7 @@ pix(tuple **     const tuples,
 int
 main(int argc, const char ** const argv) {
 
-    struct cmdlineInfo cmdline;
+    struct CmdlineInfo cmdline;
     FILE * ifP;
     struct pam inpam, outpam;
     tuple ** inTuples;
@@ -1405,8 +1412,6 @@ main(int argc, const char ** const argv) {
 
     setGlobalCP(cmdline);
 
-    srand(cmdline.randseedSpec ? cmdline.randseed : pm_randseed());
-
     ifP = pm_openr(cmdline.fileName);
 
     inTuples = pnm_readpam(ifP, &inpam, PAM_STRUCT_SIZE(tuple_type));
@@ -1421,7 +1426,8 @@ main(int argc, const char ** const argv) {
     makeAllWhite(&outpam, outTuples);
 
     if (cmdline.tri)
-        prepTrig(inpam.width, inpam.height);
+        prepTrig(inpam.width, inpam.height,
+                 cmdline.randomseedSpec, cmdline.randomseed);
     if (cmdline.quad)
         prepQuad();
 
diff --git a/editor/pbmreduce.c b/editor/pbmreduce.c
index 3a0968fe..70caa581 100644
--- a/editor/pbmreduce.c
+++ b/editor/pbmreduce.c
@@ -11,9 +11,11 @@
 */
 
 #include "pm_c_util.h"
-#include "pbm.h"
 #include "mallocvar.h"
+#include "rand.h"
 #include "shhopt.h"
+#include "pbm.h"
+
 #include <assert.h>
 
 #define SCALE 1024
@@ -105,7 +107,7 @@ parseCommandLine(int argc, const char ** argv,
         unsigned int scale;
 
         scale = strtol(argv[1], &endptr, 10);
-        if (*argv[1] == '\0') 
+        if (*argv[1] == '\0')
             pm_error("Scale argument is a null string.  Must be a number.");
         else if (*endptr != '\0')
             pm_error("Scale argument contains non-numeric character '%c'.",
@@ -115,7 +117,7 @@ parseCommandLine(int argc, const char ** argv,
                      "You specified %d", scale);
         else if (scale > INT_MAX / scale)
             pm_error("Scale argument too large.  You specified %d", scale);
-        else 
+        else
             cmdlineP->scale = scale;
 
         if (argc-1 > 1) {
@@ -145,10 +147,11 @@ struct FS {
 static void
 initializeFloydSteinberg(struct FS  * const fsP,
                          int          const newcols,
-                         unsigned int const seed,
-                         bool         const seedSpec) {
+                         bool         const seedSpec,
+                         unsigned int const seed) {
 
     unsigned int col;
+    struct pm_randSt randSt;
 
     MALLOCARRAY(fsP->thiserr, newcols + 2);
     MALLOCARRAY(fsP->nexterr, newcols + 2);
@@ -156,33 +159,36 @@ initializeFloydSteinberg(struct FS  * const fsP,
     if (fsP->thiserr == NULL || fsP->nexterr == NULL)
         pm_error("out of memory");
 
-    srand(seedSpec ? seed : pm_randseed());
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, seedSpec, seed);
 
     for (col = 0; col < newcols + 2; ++col)
-        fsP->thiserr[col] = (rand() % SCALE - HALFSCALE) / 4;
+        fsP->thiserr[col] = ((int) (pm_rand(&randSt) % SCALE) - HALFSCALE) / 4;
         /* (random errors in [-SCALE/8 .. SCALE/8]) */
+
+    pm_randterm(&randSt);
 }
 
 
 
 /*
     Scanning method
-    
+
     In Floyd-Steinberg dithering mode horizontal direction of scan alternates
     between rows; this is called "serpentine scanning".
-    
+
     Example input (14 x 7), N=3:
-    
+
     111222333444xx    Fractional pixels on the right edge and bottom edge (x)
-    111222333444xx    are ignored; their values do not influence output. 
+    111222333444xx    are ignored; their values do not influence output.
     111222333444xx
     888777666555xx
     888777666555xx
     888777666555xx
     xxxxxxxxxxxxxx
-    
+
     Output (4 x 2):
-    
+
     1234
     8765
 
@@ -196,11 +202,14 @@ enum Direction { RIGHT_TO_LEFT, LEFT_TO_RIGHT };
 static enum Direction
 oppositeDir(enum Direction const arg) {
 
+    enum Direction retval;
+
     switch (arg) {
-    case LEFT_TO_RIGHT: return RIGHT_TO_LEFT;
-    case RIGHT_TO_LEFT: return LEFT_TO_RIGHT;
+    case LEFT_TO_RIGHT: retval = RIGHT_TO_LEFT; break;
+    case RIGHT_TO_LEFT: retval = LEFT_TO_RIGHT; break;
     }
-    assert(false);  /* All cases handled above */
+
+    return retval;
 }
 
 
@@ -240,7 +249,7 @@ main(int argc, const char * argv[]) {
 
     if (cmdline.halftone == QT_FS)
         initializeFloydSteinberg(&fs, newcols,
-                                 cmdline.randomseed, cmdline.randomseedSpec);
+                                 cmdline.randomseedSpec, cmdline.randomseed);
     else {
         /* These variables are meaningless in this case, and the values
            should never be used.
@@ -258,10 +267,10 @@ main(int argc, const char * argv[]) {
         int limitCol;
         int startCol;
         int step;
-   
+
         for (colChar = 0; colChar < colChars; ++colChar)
             newbitrow[colChar] = 0x00;  /* Clear to white */
- 
+
         for (subrow = 0; subrow < cmdline.scale; ++subrow)
             pbm_readpbmrow(ifP, bitslice[subrow], cols, format);
 
@@ -274,7 +283,7 @@ main(int argc, const char * argv[]) {
         case LEFT_TO_RIGHT: {
             startCol = 0;
             limitCol = newcols;
-            step = +1;  
+            step = +1;
         } break;
         case RIGHT_TO_LEFT: {
             startCol = newcols - 1;
diff --git a/editor/pnmremap.c b/editor/pnmremap.c
index 0c0096ba..e5b59d04 100644
--- a/editor/pnmremap.c
+++ b/editor/pnmremap.c
@@ -28,6 +28,7 @@
 #include "pm_c_util.h"
 #include "mallocvar.h"
 #include "nstring.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "pam.h"
 #include "ppm.h"
@@ -399,20 +400,21 @@ randomizeError(long **       const err,
    Set a random error in the range [-1 .. 1] (normalized via FS_SCALE)
    in the error array err[][].
 -----------------------------------------------------------------------------*/
-    unsigned int const seed = (random.init == RANDOM_WITHSEED) ?
-        random.seed : pm_randseed();
-
     unsigned int col;
+    struct pm_randSt randSt;
 
     assert(random.init != RANDOM_NONE);
 
-    srand(seed);
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, random.init == RANDOM_WITHSEED, random.seed);
 
     for (col = 0; col < width; ++col) {
         unsigned int plane;
         for (plane = 0; plane < depth; ++plane)
-            err[plane][col] = rand() % (FS_SCALE * 2) - FS_SCALE;
+            err[plane][col] = pm_rand(&randSt) % (FS_SCALE * 2) - FS_SCALE;
     }
+
+    pm_randterm(&randSt);
 }
 
 
diff --git a/editor/specialty/pampaintspill.c b/editor/specialty/pampaintspill.c
index c7994723..5cd482d5 100644
--- a/editor/specialty/pampaintspill.c
+++ b/editor/specialty/pampaintspill.c
@@ -45,11 +45,11 @@
 
 #include "mallocvar.h"
 #include "nstring.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "pam.h"
 #include "pammap.h"
 
-
 static time_t const timeUpdateDelta = 30;
     /* Seconds between progress updates */
 static int const    minUpdates = 4;
@@ -67,6 +67,8 @@ struct cmdlineInfo {
     unsigned int all;
     float        power;
     unsigned int downsample;
+    unsigned int randomseedSpec;
+    unsigned int randomseed;
 };
 
 struct coords {
@@ -98,16 +100,18 @@ parseCommandLine(int argc, const char ** const argv,
     MALLOCARRAY_NOFAIL(option_def, 100);
     option_def_index = 0;          /* Incremented by OPTENTRY */
 
-    OPTENT3(0, "bgcolor",    OPT_STRING, &cmdlineP->bgcolor,    
+    OPTENT3(0, "bgcolor",    OPT_STRING, &cmdlineP->bgcolor,
             &bgcolorSpec, 0);
     OPTENT3(0, "wrap",       OPT_FLAG,   NULL,
             &cmdlineP->wrap,       0);
     OPTENT3(0, "all",        OPT_FLAG,   NULL,
             &cmdlineP->all,        0);
-    OPTENT3(0, "power",      OPT_FLOAT,  &cmdlineP->power,      
+    OPTENT3(0, "power",      OPT_FLOAT,  &cmdlineP->power,
             &powerSpec, 0);
-    OPTENT3(0, "downsample", OPT_UINT,   &cmdlineP->downsample, 
+    OPTENT3(0, "downsample", OPT_UINT,   &cmdlineP->downsample,
             &downsampleSpec, 0);
+    OPTENT3(0, "randomseed", OPT_UINT,   &cmdlineP->randomseed,
+            &cmdlineP->randomseedSpec, 0);
 
     opt.opt_table = option_def;
     opt.short_allowed = 0;
@@ -223,7 +227,9 @@ locatePaintSources(struct pam *            const pamP,
                    tuple **                const tuples,
                    tuple                   const bgColor,
                    unsigned int            const downsample,
-                   struct paintSourceSet * const paintSourcesP) {
+                   struct paintSourceSet * const paintSourcesP,
+                   bool                    const randomseedSpec,
+                   unsigned int            const randomseed) {
 /*--------------------------------------------------------------------
   Construct a list of all pixel coordinates in the input image that
   represent a non-background color.
@@ -248,21 +254,24 @@ locatePaintSources(struct pam *            const pamP,
     pm_message("Image contains %u background + %u non-background pixels",
                pamP->width * pamP->height - paintSources.size,
                paintSources.size);
-    
+
     /* Reduce the number of paint sources to reduce execution time. */
     if (downsample > 0 && downsample < paintSources.size) {
+        struct pm_randSt randSt;
         unsigned int i;
 
-        srand(pm_randseed());
+        pm_randinit(&randSt);
+        pm_srand2(&randSt, randomseedSpec, randomseed);
 
         for (i = 0; i < downsample; ++i) {
             unsigned int const swapIdx =
-                i + rand() % (paintSources.size - i);
+                i + pm_rand(&randSt) % (paintSources.size - i);
             struct coords const swapVal = paintSources.list[i];
 
             paintSources.list[i] = paintSources.list[swapIdx];
             paintSources.list[swapIdx] = swapVal;
         }
+        pm_randterm(&randSt);
         paintSources.size = downsample;
     }
 
@@ -426,7 +435,7 @@ produceOutputImage(struct pam *          const pamP,
     for (row = 0; row < pamP->height; ++row) {
         struct coords   target;
         double        * newColor;
-        
+
         MALLOCARRAY(newColor, pamP->depth);
 
         target.y = row;
@@ -484,7 +493,8 @@ main(int argc, const char *argv[]) {
                pnm_colorname(&inPam, bgColor, PAM_COLORNAME_HEXOK));
 
     locatePaintSources(&inPam, inTuples, bgColor, cmdline.downsample,
-                       &paintSources);
+                       &paintSources,
+                       cmdline.randomseedSpec, cmdline.randomseed);
 
     produceOutputImage(&inPam, inTuples, bgColor, paintSources, distFunc,
                        cmdline.power, cmdline.all, &outTuples);
@@ -498,3 +508,5 @@ main(int argc, const char *argv[]) {
 
     return 0;
 }
+
+
diff --git a/editor/specialty/ppmshift.c b/editor/specialty/ppmshift.c
index cdb0f173..85b33f3a 100644
--- a/editor/specialty/ppmshift.c
+++ b/editor/specialty/ppmshift.c
@@ -12,11 +12,11 @@
 #include <stdbool.h>
 
 #include "mallocvar.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "ppm.h"
 
 
-
 struct CmdlineInfo {
     /* All the information the user supplied in the command line,
        in a form easy for the program to use.
@@ -24,6 +24,7 @@ struct CmdlineInfo {
     const char * inputFileName;
 
     unsigned int shift;
+    unsigned int seedSpec;
     unsigned int seed;
 };
 
@@ -42,8 +43,6 @@ parseCommandLine(int argc, const char ** const argv,
 
     unsigned int option_def_index;
 
-    unsigned int seedSpec;
-
     MALLOCARRAY(option_def, 100);
 
     opt.opt_table = option_def;
@@ -52,14 +51,11 @@ parseCommandLine(int argc, const char ** const argv,
 
     option_def_index = 0;   /* incremented by OPTENT3 */
     OPTENT3(0,   "seed",            OPT_UINT,     &cmdlineP->seed,
-            &seedSpec,         0);
+            &cmdlineP->seedSpec,         0);
 
     pm_optParseOptions3(&argc, (char **)argv, opt, sizeof(opt), 0);
         /* Uses and sets argc, argv, and some of *cmdlineP and others. */
 
-    if (!seedSpec)
-        cmdlineP->seed = pm_randseed();
-
     if (argc-1 < 1)
         pm_error("You must specify the shift factor as an argument");
     else {
@@ -84,6 +80,62 @@ parseCommandLine(int argc, const char ** const argv,
 
 
 
+static void
+shiftRow(pixel *            const srcrow,
+         unsigned int       const cols,
+         unsigned int       const shift,
+         pixel *            const destrow,
+         struct pm_randSt * const randStP) {
+
+    /* the range by which a line is shifted lays in the range from */
+    /* -shift/2 .. +shift/2 pixels; however, within this range it is */
+    /* randomly chosen */
+
+    pixel * pP;
+    pixel * pP2;
+    unsigned int nowshift;
+
+    if (shift != 0)
+        nowshift = (pm_rand(randStP) % (shift+1)) - ((shift+1) / 2);
+    else
+        nowshift = 0;
+
+    pP  = &srcrow[0];
+    pP2 = &destrow[0];
+
+    /* if the shift value is less than zero, we take the original
+       pixel line and copy it into the destination line translated
+       to the left by x pixels. The empty pixels on the right end
+       of the destination line are filled up with the pixel that
+       is the right-most in the original pixel line.
+    */
+    if (nowshift < 0) {
+        unsigned int col;
+        pP += abs(nowshift);
+        for (col = 0; col < cols; ++col) {
+            PPM_ASSIGN(*pP2, PPM_GETR(*pP), PPM_GETG(*pP), PPM_GETB(*pP));
+            ++pP2;
+            if (col < (cols + nowshift) - 1)
+                ++pP;
+        }
+    } else {
+        unsigned int col;
+        /* The shift value is 0 or positive, so fill the first
+           <nowshift> pixels of the destination line with the
+           first pixel from the source line, and copy the rest of
+           the source line to the dest line
+        */
+        for (col = 0; col < cols; ++col) {
+            PPM_ASSIGN(*pP2, PPM_GETR(*pP), PPM_GETG(*pP), PPM_GETB(*pP));
+            ++pP2;
+            if (col >= nowshift)
+                ++pP;
+        }
+    }
+}
+
+
+
 int
 main(int argc, const char ** argv) {
 
@@ -95,13 +147,15 @@ main(int argc, const char ** argv) {
     pixel * destrow;
     unsigned int row;
     unsigned int shift;
+    struct pm_randSt randSt;
 
     /* parse in 'default' parameters */
     pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
 
-    srand(cmdline.seed);
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, cmdline.seedSpec, cmdline.seed);
 
     ifP = pm_openr(cmdline.inputFileName);
 
@@ -115,67 +169,23 @@ main(int argc, const char ** argv) {
     } else
         shift = cmdline.shift;
 
-    srcrow = ppm_allocrow(cols);
-
+    srcrow  = ppm_allocrow(cols);
     destrow = ppm_allocrow(cols);
 
     ppm_writeppminit(stdout, cols, rows, maxval, 0);
 
-    /** now do the shifting **/
-    /* the range by which a line is shifted lays in the range from */
-    /* -shift/2 .. +shift/2 pixels; however, within this range it is */
-    /* randomly chosen */
     for (row = 0; row < rows; ++row) {
-        pixel * pP;
-        pixel * pP2;
-        unsigned int nowshift;
-
-        if (shift != 0)
-            nowshift = (rand() % (shift+1)) - ((shift+1) / 2);
-        else
-            nowshift = 0;
-
         ppm_readppmrow(ifP, srcrow, cols, maxval, format);
 
-        pP  = &srcrow[0];
-        pP2 = &destrow[0];
-
-        /* if the shift value is less than zero, we take the original
-           pixel line and copy it into the destination line translated
-           to the left by x pixels. The empty pixels on the right end
-           of the destination line are filled up with the pixel that
-           is the right-most in the original pixel line.
-        */
-        if (nowshift < 0) {
-            unsigned int col;
-            pP += abs(nowshift);
-            for (col = 0; col < cols; ++col) {
-                PPM_ASSIGN(*pP2, PPM_GETR(*pP), PPM_GETG(*pP), PPM_GETB(*pP));
-                ++pP2;
-                if (col < (cols + nowshift) - 1)
-                    ++pP;
-            }
-        } else {
-            unsigned int col;
-            /* The shift value is 0 or positive, so fill the first
-               <nowshift> pixels of the destination line with the
-               first pixel from the source line, and copy the rest of
-               the source line to the dest line
-            */
-            for (col = 0; col < cols; ++col) {
-                PPM_ASSIGN(*pP2, PPM_GETR(*pP), PPM_GETG(*pP), PPM_GETB(*pP));
-                ++pP2;
-                if (col >= nowshift)
-                    ++pP;
-            }
-        }
+        shiftRow(srcrow, cols, shift, destrow, &randSt);
 
         ppm_writeppmrow(stdout, destrow, cols, maxval, 0);
     }
 
-    pm_close(ifP);
-    ppm_freerow(srcrow);
     ppm_freerow(destrow);
+    ppm_freerow(srcrow);
+    pm_close(ifP);
+    pm_randterm(&randSt);
 
     return 0;
 }
diff --git a/editor/specialty/ppmspread.c b/editor/specialty/ppmspread.c
index 6753f4fe..7b9558e3 100644
--- a/editor/specialty/ppmspread.c
+++ b/editor/specialty/ppmspread.c
@@ -10,102 +10,158 @@
 
 #include <string.h>
 
+#include "nstring.h"
+#include "rand.h"
+#include "shhopt.h"
 #include "ppm.h"
 
 
+struct CmdlineInfo {
+    /* This structure represents all of the information the user
+       supplied in the command line but in a form that's easy for the
+       program to use.
+    */
+    const char * inputFilename;  /* '-' if stdin */
+    unsigned int spread;
+    unsigned int randomseedSpec;
+    unsigned int randomseed;
+};
+
+
+
+static void
+parseCommandLine(int argc, const char ** const argv,
+                 struct CmdlineInfo * const cmdlineP ) {
+
+    optEntry     * option_def;
+        /* Instructions to OptParseOptions3 on how to parse our options */
+    optStruct3     opt;
+    unsigned int   option_def_index;
+
+    MALLOCARRAY_NOFAIL(option_def, 100);
+    option_def_index = 0;          /* Incremented by OPTENTRY */
+
+    OPTENT3(0, "randomseed", OPT_UINT,   &cmdlineP->randomseed,
+            &cmdlineP->randomseedSpec, 0);
+
+    opt.opt_table = option_def;
+    opt.short_allowed = 0;
+    opt.allowNegNum = 1;
+
+    pm_optParseOptions3( &argc, (char **)argv, opt, sizeof(opt), 0 );
+
+    if (argc-1 < 1)
+        pm_error("You must specify the spread factor as an argument");
+    else {
+        const char * error;
+        pm_string_to_uint(argv[1], &cmdlineP->spread, &error);
+
+        if (error)
+            pm_error("Spread factor '%s' is not an unsigned integer.  %s",
+                     argv[1], error);
+
+        if (argc-1 < 2)
+            cmdlineP->inputFilename = "-";
+        else {
+            cmdlineP->inputFilename = argv[2];
+            if (argc-1 >2)
+                pm_error("Too many arguments: %u.  "
+                         "The only possible arguments are "
+                         "the spread factor and the optional input file name",
+                         argc-1);
+        }
+    }
+}
+
+
+
+static void
+spreadRow(pixel **           const srcarray,
+          unsigned int       const cols,
+          unsigned int       const rows,
+          unsigned int       const spread,
+          unsigned int       const row,
+          pixel **           const destarray,
+          struct pm_randSt * const randStP) {
+
+    unsigned int col;
+
+    for (col = 0; col < cols; ++col) {
+        pixel const p = srcarray[row][col];
+
+        int const xdis = (pm_rand(randStP) % (spread + 1) )
+            - ((spread + 1) / 2);
+        int const ydis = (pm_rand(randStP) % (spread + 1))
+            - ((spread + 1) / 2);
+
+        int const xnew = col + xdis;
+        int const ynew = row + ydis;
+
+        /* only set the displaced pixel if it's within the bounds
+           of the image
+        */
+        if (xnew >= 0 && xnew < cols && ynew >= 0 && ynew < rows) {
+            /* Displacing a pixel is accomplished by swapping it
+               with another pixel in its vicinity.
+            */
+            pixel const p2 = srcarray[ynew][xnew];
+                /* Original value of second pixel */
+
+            /* Set second pixel to new value */
+            PPM_ASSIGN(destarray[ynew][xnew],
+                       PPM_GETR(p), PPM_GETG(p), PPM_GETB(p));
+
+            /* Set first pixel to (old) value of second */
+            PPM_ASSIGN(destarray[row][col],
+                       PPM_GETR(p2), PPM_GETG(p2), PPM_GETB(p2));
+        } else {
+            /* Displaced pixel is out of bounds; leave the old pixel there.
+            */
+            PPM_ASSIGN(destarray[row][col],
+                       PPM_GETR(p), PPM_GETG(p), PPM_GETB(p));
+        }
+    }
+}
+
+
 
 int
-main(int    argc,
-     char * argv[]) {
+main(int          argc,
+     const char * argv[]) {
 
+    struct CmdlineInfo cmdline;
     FILE * ifP;
-    int argn, rows, cols;
+    int rows, cols;
     unsigned int row;
-    pixel ** destarray, ** srcarray;
-    pixel * pP;
-    pixel * pP2;
+    pixel ** destarray;
+    pixel ** srcarray;
     pixval maxval;
-    pixval r1, g1, b1;
-    int amount;
-    const char * const usage = "amount [ppmfile]\n        amount: # of pixels to displace a pixel by at most\n";
-
-    /* parse in 'default' parameters */
-    ppm_init(&argc, argv);
-
-    argn = 1;
-
-    /* parse in amount & seed */
-    if (argn == argc)
-        pm_usage(usage);
-    if (sscanf(argv[argn], "%d", &amount) != 1)
-        pm_usage(usage);
-    if (amount < 0)
-        pm_error("amount should be a positive number");
-    ++argn;
-
-    /* parse in filename (if present, stdin otherwise) */
-    if (argn != argc)
-    {
-        ifP = pm_openr(argv[argn]);
-        ++argn;
-    }
-    else
-        ifP = stdin;
+    struct pm_randSt randSt;
 
-    if (argn != argc)
-        pm_usage(usage);
+    pm_proginit(&argc, argv);
+
+    parseCommandLine(argc, argv, &cmdline);
+
+    ifP = pm_openr(cmdline.inputFilename);
 
     srcarray = ppm_readppm(ifP, &cols, &rows, &maxval);
 
     destarray = ppm_allocarray(cols, rows);
 
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, cmdline.randomseedSpec, cmdline.randomseed);
+
     /* clear out the buffer */
     for (row = 0; row < rows; ++row)
         memset(destarray[row], 0, cols * sizeof(pixel));
 
-    srand(pm_randseed());
-
-    /* start displacing pixels */
+    /* Displace pixels */
     for (row = 0; row < rows; ++row) {
-        unsigned int col;
-        pP = &srcarray[row][0];
-
-        for (col = 0; col < cols; ++col) {
-            int const xdis = (rand() % (amount+1)) - ((amount+1) / 2);
-            int const ydis = (rand() % (amount+1)) - ((amount+1) / 2);
+        spreadRow(srcarray, cols, rows, cmdline.spread, row,
+                  destarray, &randSt);
 
-            int const xnew = col + xdis;
-            int const ynew = row + ydis;
-
-            /* only set the displaced pixel if it's within the bounds
-               of the image
-            */
-            if (xnew >= 0 && xnew < cols && ynew >= 0 && ynew < rows) {
-                /* displacing a pixel is accomplished by swapping it
-                   with another pixel in its vicinity - so, first
-                   store other pixel's RGB
-                */
-                pP2 = &srcarray[ynew][xnew];
-                r1 = PPM_GETR(*pP2);
-                g1 = PPM_GETG(*pP2);
-                b1 = PPM_GETB(*pP2);
-                /* set second pixel to new value */
-                pP2 = &destarray[ynew][xnew];
-                PPM_ASSIGN(*pP2, PPM_GETR(*pP), PPM_GETG(*pP), PPM_GETB(*pP));
-                
-                /* now, set first pixel to (old) value of second */
-                pP2 = &destarray[row][col];
-                PPM_ASSIGN(*pP2, r1, g1, b1);
-            } else {
-                /* displaced pixel is out of bounds; leave the old
-                   pixel there
-                */
-                pP2 = &destarray[row][col];
-                PPM_ASSIGN(*pP2, PPM_GETR(*pP), PPM_GETG(*pP), PPM_GETB(*pP));
-            }
-            ++pP;
-        }
     }
+    pm_randterm(&randSt);
 
     ppm_writeppm(stdout, destarray, cols, rows, maxval, 0);
 
@@ -115,3 +171,5 @@ main(int    argc,
 
     return 0;
 }
+
+
diff --git a/generator/pamcrater.c b/generator/pamcrater.c
index 43c27dbc..b8ceafa5 100644
--- a/generator/pamcrater.c
+++ b/generator/pamcrater.c
@@ -48,6 +48,7 @@
 
 #include "pm_c_util.h"
 #include "mallocvar.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "nstring.h"
 #include "pam.h"
@@ -169,11 +170,12 @@ static double const DepthBias2  = 0.5;      /* Square of depth bias */
 
 
 static double const
-cast(double const high) {
+cast(double             const high,
+     struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    A random number in the range [0, 'high'].
 -----------------------------------------------------------------------------*/
-    return high * ((rand() & 0x7FFF) / arand);
+    return high * ((pm_rand(randStP) & 0x7FFF) / arand);
 }
 
 
@@ -252,11 +254,12 @@ setElev(struct pam * const pamP,
 
 
 static void
-smallCrater(struct pam * const pamP,
-            tuple **     const terrain,
-            int          const cx,
-            int          const cy,
-            double       const radius) {
+smallCrater(struct pam *       const pamP,
+            tuple **           const terrain,
+            int                const cx,
+            int                const cy,
+            double             const radius,
+            struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
    Generate a crater with a special method for tiny craters.
 
@@ -283,10 +286,10 @@ smallCrater(struct pam * const pamP,
             /* The mean elevation of the Moore neighborhood (9 pixels
                centered on the crater location).
             */
-        
+
         /* Perturb the mean elevation by a small random factor. */
 
-        int const x = radius >= 1 ? ((rand() >> 8) & 0x3) - 1 : 0;
+        int const x = radius >= 1 ? ((pm_rand(randStP) >> 8) & 0x3) - 1 : 0;
 
         assert(axelev > 0);
 
@@ -374,7 +377,7 @@ normalCrater(struct pam * const pamP,
                 av = (axelev + cz) * (1 - roll) +
                     (terrainMod(pamP, terrain, x, y) + cz) * roll;
                 av = MAX(1000, MIN(64000, av));
-                
+
                 setElev(pamP, terrain, x, y, av);
             }
         }
@@ -388,19 +391,20 @@ normalCrater(struct pam * const pamP,
 
 
 static void
-plopCrater(struct pam * const pamP,
-           tuple **     const terrain,
-           int          const cx,
-           int          const cy,
-           double       const radius,
-           bool         const verbose) {
+plopCrater(struct pam *       const pamP,
+           tuple **           const terrain,
+           int                const cx,
+           int                const cy,
+           double             const radius,
+           bool               const verbose,
+           struct pm_randSt * const randStP) {
 
     if (verbose && pm_have_float_format())
         pm_message("Plopping crater at (%4d, %4d) with radius %g",
                    cx, cy, radius);
 
     if (radius < 3)
-        smallCrater (pamP, terrain, cx, cy, radius);
+        smallCrater (pamP, terrain, cx, cy, radius, randStP);
     else
         normalCrater(pamP, terrain, cx, cy, radius);
 }
@@ -448,6 +452,7 @@ genCraters(struct CmdlineInfo const cmdline) {
 -----------------------------------------------------------------------------*/
     tuple ** terrain;    /* elevation array */
     struct pam pam;
+    struct pm_randSt randSt;
 
     /* Allocate the elevation array and initialize it to mean surface
        elevation.
@@ -455,18 +460,22 @@ genCraters(struct CmdlineInfo const cmdline) {
 
     initCanvas(cmdline.width, cmdline.height, &pam, &terrain);
 
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, cmdline.randomseedSpec, cmdline.randomseed);
+
     if (cmdline.test)
         plopCrater(&pam, terrain,
                    pam.width/2 + cmdline.offset,
                    pam.height/2 + cmdline.offset,
-                   (double) cmdline.radius, cmdline.verbose);
+                   (double) cmdline.radius, cmdline.verbose,
+                   &randSt);
     else {
         unsigned int const ncraters = cmdline.number; /* num of craters */
         unsigned int l;
 
         for (l = 0; l < ncraters; ++l) {
-            int const cx = cast((double) pam.width  - 1);
-            int const cy = cast((double) pam.height - 1);
+            int const cx = cast((double) pam.width  - 1, &randSt);
+            int const cy = cast((double) pam.height - 1, &randSt);
 
             /* Thanks, Rudy, for this equation that maps the uniformly
                distributed numbers from cast() into an area-law distribution
@@ -475,9 +484,11 @@ genCraters(struct CmdlineInfo const cmdline) {
                Produces values within the interval:
                0.56419 <= radius <= 56.419
             */
-            double const radius = sqrt(1 / (M_PI * (1 - cast(0.9999))));
+            double const radius =
+                sqrt(1 / (M_PI * (1 - cast(0.9999, &randSt))));
 
-            plopCrater(&pam, terrain, cx, cy, radius, cmdline.verbose);
+            plopCrater(&pam, terrain, cx, cy, radius,
+                       cmdline.verbose, &randSt);
 
             if (((l + 1) % 100000) == 0)
                 pm_message("%u craters generated of %u (%u%% done)",
@@ -485,6 +496,8 @@ genCraters(struct CmdlineInfo const cmdline) {
         }
     }
 
+    pm_randterm(&randSt);
+
     pnm_writepam(&pam, terrain);
 
     pnm_freepamarray(terrain, &pam);
@@ -503,12 +516,9 @@ main(int argc, const char ** argv) {
 
     parseCommandLine(argc, argv, &cmdline);
 
-    srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
-
     genCraters(cmdline);
 
     return 0;
 }
 
 
-
diff --git a/generator/pamstereogram.c b/generator/pamstereogram.c
index 39356294..6e64b127 100644
--- a/generator/pamstereogram.c
+++ b/generator/pamstereogram.c
@@ -58,11 +58,12 @@
 #include "pm_c_util.h"
 #include "mallocvar.h"
 #include "nstring.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "pam.h"
 
 
-enum outputType {OUTPUT_BW, OUTPUT_GRAYSCALE, OUTPUT_COLOR};
+enum OutputType {OUTPUT_BW, OUTPUT_GRAYSCALE, OUTPUT_COLOR};
 
 /* ---------------------------------------------------------------------- */
 
@@ -92,7 +93,7 @@ struct cmdlineInfo {
     unsigned int smoothing;      /* -smoothing option */
     unsigned int randomseed;     /* -randomseed option */
     unsigned int randomseedSpec; /* -randomseed option count */
-    enum outputType outputType;  /* Type of output file */
+    enum OutputType outputType;  /* Type of output file */
     unsigned int xbegin;         /* -xbegin option */
     unsigned int xbeginSpec;     /* -xbegin option count */
     unsigned int tileable;       /* -tileable option */
@@ -406,11 +407,12 @@ typedef struct outGenerator {
 
 
 
-struct randomState {
+struct RandomState {
     /* The state of a randomColor generator. */
     unsigned int magnifypat;
     tuple *      currentRow;
     unsigned int prevy;
+    struct pm_randSt * randStP;
 };
 
 
@@ -425,7 +427,7 @@ randomColor(outGenerator * const outGenP,
 /*----------------------------------------------------------------------------
    Return a random RGB value.
 -----------------------------------------------------------------------------*/
-    struct randomState * const stateP = outGenP->stateP;
+    struct RandomState * const stateP = outGenP->stateP;
 
     /* Every time we start a new row, we select a new sequence of random
        colors.
@@ -440,7 +442,7 @@ randomColor(outGenerator * const outGenP,
             unsigned int plane;
 
             for (plane = 0; plane < outGenP->pam.depth; ++plane) {
-                unsigned int const randval = rand();
+                unsigned int const randval = pm_rand(stateP->randStP);
                 thisTuple[plane] = randval % modulus;
             }
         }
@@ -460,9 +462,13 @@ static outGenStateTerm termRandomColor;
 static void
 termRandomColor(outGenerator * const outGenP) {
 
-    struct randomState * const stateP = outGenP->stateP;
+    struct RandomState * const stateP = outGenP->stateP;
 
     pnm_freepamrow(stateP->currentRow);
+
+    pm_randterm(stateP->randStP);
+
+    free(stateP->randStP);
 }
 
 
@@ -472,7 +478,7 @@ initRandomColor(outGenerator *     const outGenP,
                 const struct pam * const inPamP,
                 struct cmdlineInfo const cmdline) {
 
-    struct randomState * stateP;
+    struct RandomState * stateP;
 
     outGenP->pam.format      = PAM_FORMAT;
     outGenP->pam.plainformat = 0;
@@ -503,6 +509,10 @@ initRandomColor(outGenerator *     const outGenP,
     stateP->magnifypat = cmdline.magnifypat;
     stateP->prevy      = (unsigned int)(-cmdline.magnifypat);
 
+    MALLOCVAR_NOFAIL(stateP->randStP);
+    pm_randinit(stateP->randStP);
+    pm_srand2(stateP->randStP, cmdline.randomseedSpec, cmdline.randomseed);
+
     outGenP->stateP         = stateP;
     outGenP->getTuple       = &randomColor;
     outGenP->terminateState = &termRandomColor;
@@ -1600,8 +1610,6 @@ main(int argc, const char *argv[]) {
     if (cmdline.verbose)
         reportParameters(cmdline);
 
-    srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
-
     ifP = pm_openr(cmdline.inputFilespec);
 
     /* Produce a stereogram. */
@@ -1613,4 +1621,3 @@ main(int argc, const char *argv[]) {
 }
 
 
-
diff --git a/generator/pgmnoise.c b/generator/pgmnoise.c
index e3e4eca6..f5f48ff3 100644
--- a/generator/pgmnoise.c
+++ b/generator/pgmnoise.c
@@ -7,12 +7,12 @@
 #include "pm_c_util.h"
 #include "mallocvar.h"
 #include "nstring.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "pgm.h"
 
 
-
-struct cmdlineInfo {
+struct CmdlineInfo {
     /* All the information the user supplied in the command line,
        in a form easy for the program to use.
     */
@@ -21,13 +21,15 @@ struct cmdlineInfo {
     unsigned int maxval;
     unsigned int randomseed;
     unsigned int randomseedSpec;
+    unsigned int verbose;
 };
 
 
 
 static void
-parseCommandLine(int argc, const char ** const argv,
-                 struct cmdlineInfo * const cmdlineP) {
+parseCommandLine(int argc,
+                 const char **        const argv,
+                 struct CmdlineInfo * const cmdlineP) {
 /*----------------------------------------------------------------------------
    Note that the file spec array we return is stored in the storage that
    was passed to us as the argv array.
@@ -46,6 +48,8 @@ parseCommandLine(int argc, const char ** const argv,
             &cmdlineP->randomseedSpec,      0);
     OPTENT3(0,   "maxval",       OPT_UINT,    &cmdlineP->maxval,
             &maxvalSpec,                    0);
+    OPTENT3(0,   "verbose",      OPT_FLAG,    NULL,
+            &cmdlineP->verbose,             0);
 
     opt.opt_table = option_def;
     opt.short_allowed = FALSE;  /* We have no short (old-fashioned) options */
@@ -57,7 +61,7 @@ parseCommandLine(int argc, const char ** const argv,
 
     if (maxvalSpec) {
         if (cmdlineP->maxval > PGM_OVERALLMAXVAL)
-            pm_error("Maxval too large: %u.  Maximu is %u",
+            pm_error("Maxval too large: %u.  Maximum is %u",
                      cmdlineP->maxval, PGM_OVERALLMAXVAL);
         else if (cmdlineP->maxval == 0)
             pm_error("Maxval must not be zero");
@@ -95,15 +99,21 @@ parseCommandLine(int argc, const char ** const argv,
 
 
 static unsigned int
-randPool(unsigned int const digits) {
+randPool(unsigned int       const digits,
+         bool               const verbose,
+         struct pm_randSt * const randStP) {
 /*----------------------------------------------------------------------------
   Draw 'digits' bits from pool of random bits.  If the number of random bits
-  in pool is insufficient, call rand() and add 31 bits to it.
+  in pool is insufficient, call pm_rand() and add N bits to it.
+
+  N is 31 or 32.  In raw mode we set this value to 32 regardless of the
+  actual number of available bits.  If insufficient, the MSB will be
+  set to zero.
 
   'digits' must be at most 16.
 
-  We assume that each call to rand() generates 31 bits, or RAND_MAX ==
-  2147483647.
+  We assume that each call to pm_rand() generates 31 or 32 bits,
+  or randStP->max == 2147483647 or 4294967294.
 
   The underlying logic is flexible and endian-free.  The above conditions
   can be relaxed.
@@ -112,21 +122,24 @@ randPool(unsigned int const digits) {
     static unsigned int len=0;        /* number of valid bits in pool */
 
     unsigned int const mask = (1 << digits) - 1;
-
+    unsigned int const randbits = (randStP->max == 2147483647) ? 31 : 32;
     unsigned int retval;
 
-    assert(RAND_MAX == 2147483647 && digits <= 16);
+    assert(randStP->max == 2147483647 || randStP->max == 4294967294);
+    assert(digits <= 16);
 
     retval = hold;  /* initial value */
 
     if (len > digits) { /* Enough bits in hold to satisfy request */
         hold >>= digits;
         len   -= digits;
-    } else {              /* Load another 31 bits into hold */
-        hold    = rand();
+    } else {            /* Load another 31 or 32 bits into hold */
+        hold    = pm_rand(randStP);
+        if (verbose)
+            pm_message("pm_rand(): %08lX", hold);
         retval |= (hold << len);
         hold >>=  (digits - len);
-        len = 31 - digits + len;
+        len = randbits - digits + len;
     }
     return (retval & mask);
 }
@@ -134,19 +147,23 @@ randPool(unsigned int const digits) {
 
 
 static void
-pgmnoise(FILE *       const ofP,
-         unsigned int const cols,
-         unsigned int const rows,
-         gray         const maxval) {
-
-    bool const usingPool = !(RAND_MAX==2147483647 && (maxval & (maxval+1)));
+pgmnoise(FILE *             const ofP,
+         unsigned int       const cols,
+         unsigned int       const rows,
+         gray               const maxval,
+         bool               const verbose,
+         struct pm_randSt * const randStP) {
+
+    bool const usingPool =
+      !( (randStP->max==2147483647UL || randStP->max==4294967294UL)
+           && (maxval & (maxval+1)) );
     unsigned int const bitLen = pm_maxvaltobits(maxval);
 
     unsigned int row;
     gray * destrow;
 
     /* If maxval is 2^n-1, we draw exactly n bits from the pool.
-       Otherwise call rand() and determine gray value by modulo.
+       Otherwise call pm_rand() and determine gray value by modulo.
 
        In the latter case, there is a minuscule skew toward 0 (=black)
        because smaller numbers are produced more frequently by modulo.
@@ -155,7 +172,7 @@ pgmnoise(FILE *       const ofP,
 
        To illustrate the point, consider converting the outcome of one
        roll of a fair, six-sided die to 5 values (0 to 4) by N % 5.  The
-       probability for values 1, 2, 3, 4 are 1/6, but 0 alone is 2/6.
+       probability for values 1, 2, 3, 4 is 1/6, but 0 alone is 2/6.
        Average is 10/6 or 1.6667, compared to 2.0 from an ideal
        generator which produces exactly 5 values.  With two dice
        average improves to 70/36 or 1.9444.
@@ -172,12 +189,11 @@ pgmnoise(FILE *       const ofP,
         if (usingPool) {
             unsigned int col;
             for (col = 0; col < cols; ++col)
-                destrow[col] = randPool(bitLen);
-        }
-        else {
+                destrow[col] = randPool(bitLen, verbose, randStP);
+        } else {
             unsigned int col;
             for (col = 0; col < cols; ++col)
-                destrow[col] = rand() % (maxval + 1);
+                destrow[col] = pm_rand(randStP) % (maxval + 1);
         }
         pgm_writepgmrow(ofP, destrow, cols, maxval, 0);
     }
@@ -191,18 +207,22 @@ int
 main(int          argc,
      const char * argv[]) {
 
-    struct cmdlineInfo cmdline;
+    struct CmdlineInfo cmdline;
+    struct pm_randSt randSt;
 
     pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
 
-    srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, cmdline.randomseedSpec, cmdline.randomseed);
+
+    pgmnoise(stdout, cmdline.width, cmdline.height, cmdline.maxval,
+             cmdline.verbose, &randSt);
 
-    pgmnoise(stdout, cmdline.width, cmdline.height, cmdline.maxval);
+    pm_randterm(&randSt);
 
     return 0;
 }
 
 
-
diff --git a/generator/ppmforge.c b/generator/ppmforge.c
index 114f7f18..c7cfdb84 100644
--- a/generator/ppmforge.c
+++ b/generator/ppmforge.c
@@ -37,9 +37,10 @@
 #include <assert.h>
 
 #include "pm_c_util.h"
-#include "ppm.h"
 #include "mallocvar.h"
+#include "rand.h"
 #include "shhopt.h"
+#include "ppm.h"
 
 static double const hugeVal = 1e50;
 
@@ -59,12 +60,8 @@ typedef struct {
 
 /* Definition for obtaining random numbers. */
 
-#define nrand 4               /* Gauss() sample count */
-#define Cast(low, high) ((low)+(((high)-(low)) * ((rand() & 0x7FFF) / arand)))
-
 /*  Local variables  */
 
-static double arand, gaussadd, gaussfac; /* Gaussian random parameters */
 static double fracdim;            /* Fractal dimension */
 static double powscale;           /* Power law scaling exponent */
 static int meshsize = 256;        /* FFT mesh size */
@@ -335,45 +332,73 @@ fourn(float *     const data,
 #undef SWAP
 
 
+struct Gauss {
+    struct pm_randSt randSt;
+    unsigned int     nrand;          /* Gauss() sample count */
+    double           arand;
+    double           gaussadd;
+    double           gaussfac;
+};
+
+
 
 static void
-initgauss(unsigned int const seed) {
+initgauss(struct Gauss * const gaussP,
+          unsigned int   const seed) {
 /*----------------------------------------------------------------------------
   Initialize random number generators.
 
   As given in Peitgen & Saupe, page 77.
 -----------------------------------------------------------------------------*/
+    gaussP->nrand = 4;
+
     /* Range of random generator */
-    arand = pow(2.0, 15.0) - 1.0;
-    gaussadd = sqrt(3.0 * nrand);
-    gaussfac = 2 * gaussadd / (nrand * arand);
-    srand(seed);
+    gaussP->arand    = pow(2.0, 15.0) - 1.0;
+    gaussP->gaussadd = sqrt(3.0 * gaussP->nrand);
+    gaussP->gaussfac = 2 * gaussP->gaussadd / (gaussP->nrand * gaussP->arand);
+
+    pm_randinit(&gaussP->randSt);
+    pm_srand(&gaussP->randSt, seed);
 }
 
 
 
 static double
-gauss() {
+gauss(struct Gauss * const gaussP) {
 /*----------------------------------------------------------------------------
   A Gaussian random number.
 
   As given in Peitgen & Saupe, page 77.
 -----------------------------------------------------------------------------*/
-    int i;
+    unsigned int i;
     double sum;
 
-    for (i = 1, sum = 0.0; i <= nrand; ++i) {
-        sum += (rand() & 0x7FFF);
+    for (i = 1, sum = 0.0; i <= gaussP->nrand; ++i) {
+        sum += (pm_rand(&gaussP->randSt) & 0x7FFF);
     }
-    return gaussfac * sum - gaussadd;
+    return gaussP->gaussfac * sum - gaussP->gaussadd;
+}
+
+
+
+static double
+cast(double         const low,
+     double         const high,
+     struct Gauss * const gaussP) {
+
+    return
+        low +
+        ((high-low) * ((pm_rand(&gaussP->randSt) & 0x7FFF) / gaussP->arand));
+
 }
 
 
 
 static void
-spectralsynth(float **     const x,
-              unsigned int const n,
-              double        const h) {
+spectralsynth(float **       const x,
+              unsigned int   const n,
+              double         const h,
+              struct Gauss * const gaussP) {
 /*----------------------------------------------------------------------------
   Spectrally synthesized fractal motion in two dimensions.
 
@@ -395,9 +420,10 @@ spectralsynth(float **     const x,
 
     for (i = 0; i <= n / 2; i++) {
         for (j = 0; j <= n / 2; j++) {
-            phase = 2 * M_PI * ((rand() & 0x7FFF) / arand);
+            phase = 2 * M_PI * ((pm_rand(&gaussP->randSt) & 0x7FFF) /
+                                gaussP->arand);
             if (i != 0 || j != 0) {
-                rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss();
+                rad = pow((double) (i*i + j*j), -(h + 1) / 2) * gauss(gaussP);
             } else {
                 rad = 0;
             }
@@ -416,8 +442,9 @@ spectralsynth(float **     const x,
     Imag(a, n / 2, n / 2) = 0;
     for (i = 1; i <= n / 2 - 1; i++) {
         for (j = 1; j <= n / 2 - 1; j++) {
-            phase = 2 * M_PI * ((rand() & 0x7FFF) / arand);
-            rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss();
+            phase = 2 * M_PI * ((pm_rand(&gaussP->randSt) & 0x7FFF) /
+                                gaussP->arand);
+            rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss(gaussP);
             rcos = rad * cos(phase);
             rsin = rad * sin(phase);
             Real(a, i, n - j) = rcos;
@@ -469,22 +496,22 @@ temprgb(double   const temp,
 
 
 static void
-etoile(pixel * const pix) {
+etoile(pixel *        const pix,
+       struct Gauss * const gaussP) {
 /*----------------------------------------------------------------------------
     Set a pixel in the starry sky.
 -----------------------------------------------------------------------------*/
-    if ((rand() % 1000) < starfraction) {
-#define StarQuality 0.5       /* Brightness distribution exponent */
-#define StarIntensity   8         /* Brightness scale factor */
-#define StarTintExp 0.5       /* Tint distribution exponent */
-        double v = StarIntensity * pow(1 / (1 - Cast(0, 0.9999)),
-                                       (double) StarQuality),
-            temp,
-            r, g, b;
-
-        if (v > 255) {
-            v = 255;
-        }
+    if ((pm_rand(&gaussP->randSt) % 1000) < starfraction) {
+        double const starQuality   = 0.5;
+            /* Brightness distribution exponent */
+        double const starIntensity = 8;
+            /* Brightness scale factor */
+        double const starTintExp = 0.5;
+            /* Tint distribution exponent */
+        double const v =
+            MIN(255.0,
+                starIntensity * pow(1 / (1 - cast(0, 0.9999, gaussP)),
+                                    (double) starQuality));
 
         /* We make a special case for star color  of zero in order to
            prevent  floating  point  roundoff  which  would  otherwise
@@ -493,13 +520,17 @@ etoile(pixel * const pix) {
            256 shades in the image. */
 
         if (starcolor == 0) {
-            int vi = v;
+            pixval const vi = v;
 
             PPM_ASSIGN(*pix, vi, vi, vi);
         } else {
+            double temp;
+            double r, g, b;
+
             temp = 5500 + starcolor *
-                pow(1 / (1 - Cast(0, 0.9999)), StarTintExp) *
-                ((rand() & 7) ? -1 : 1);
+                pow(1 / (1 - cast(0, 0.9999, gaussP)), starTintExp) *
+                ((pm_rand(&gaussP->randSt) & 7) ? -1 : 1);
+
             /* Constrain temperature to a reasonable value: >= 2600K
                (S Cephei/R Andromedae), <= 28,000 (Spica). */
             temp = MAX(2600, MIN(28000, temp));
@@ -575,7 +606,8 @@ createPlanetStuff(bool             const clouds,
                   unsigned char ** const cpP,
                   Vector *         const sunvecP,
                   unsigned int     const cols,
-                  pixval           const maxval) {
+                  pixval           const maxval,
+                  struct Gauss *   const gaussP) {
 
     double *u, *u1;
     unsigned int *bxf, *bxc;
@@ -585,8 +617,8 @@ createPlanetStuff(bool             const clouds,
 
     /* Compute incident light direction vector. */
 
-    shang = hourspec ? hourangle : Cast(0, 2 * M_PI);
-    siang = inclspec ? inclangle : Cast(-M_PI * 0.12, M_PI * 0.12);
+    shang = hourspec ? hourangle : cast(0, 2 * M_PI, gaussP);
+    siang = inclspec ? inclangle : cast(-M_PI * 0.12, M_PI * 0.12, gaussP);
 
     sunvecP->x = sin(shang) * cos(siang);
     sunvecP->y = sin(siang);
@@ -594,7 +626,7 @@ createPlanetStuff(bool             const clouds,
 
     /* Allow only 25% of random pictures to be crescents */
 
-    if (!hourspec && ((rand() % 100) < 75)) {
+    if (!hourspec && ((pm_rand(&gaussP->randSt) % 100) < 75)) {
         flipped = (sunvecP->z < 0);
         sunvecP->z = fabs(sunvecP->z);
     } else
@@ -634,7 +666,7 @@ createPlanetStuff(bool             const clouds,
         pm_error("Cannot allocate %u element interpolation tables.", cols);
     {
         unsigned int j;
-        for (j = 0; j < cols; j++) {
+        for (j = 0; j < cols; ++j) {
             double const bx = (n - 1) * uprj(j, cols);
 
             bxf[j] = floor(bx);
@@ -651,8 +683,9 @@ createPlanetStuff(bool             const clouds,
 
 
 static void
-generateStarrySkyRow(pixel *      const pixels,
-                     unsigned int const cols) {
+generateStarrySkyRow(pixel *        const pixels,
+                     unsigned int   const cols,
+                     struct Gauss * const gaussP) {
 /*----------------------------------------------------------------------------
   Generate a starry sky.  Note  that no FFT is performed;
   the output is  generated  directly  from  a  power  law
@@ -660,8 +693,8 @@ generateStarrySkyRow(pixel *      const pixels,
 -----------------------------------------------------------------------------*/
     unsigned int j;
 
-    for (j = 0; j < cols; j++)
-        etoile(pixels + j);
+    for (j = 0; j < cols; ++j)
+        etoile(pixels + j, gaussP);
 }
 
 
@@ -881,7 +914,8 @@ generatePlanetRow(pixel *         const pixelrow,
                   int             const byc,
                   int             const byf,
                   Vector          const sunvec,
-                  pixval          const maxval) {
+                  pixval          const maxval,
+                  struct Gauss *  const gaussP) {
 
     unsigned int const StarClose = 2;
 
@@ -924,24 +958,25 @@ generatePlanetRow(pixel *         const pixelrow,
     /* Left stars */
 
     for (col = 0; (int)col < (int)(cols/2 - (lcos + StarClose)); ++col)
-        etoile(&pixelrow[col]);
+        etoile(&pixelrow[col], gaussP);
 
     /* Right stars */
 
     for (col = cols/2 + (lcos + StarClose); col < cols; ++col)
-        etoile(&pixelrow[col]);
+        etoile(&pixelrow[col], gaussP);
 }
 
 
 
 static void
-genplanet(bool         const stars,
-          bool         const clouds,
-          float *      const a,
-          unsigned int const cols,
-          unsigned int const rows,
-          unsigned int const n,
-          unsigned int const rseed) {
+genplanet(bool           const stars,
+          bool           const clouds,
+          float *        const a,
+          unsigned int   const cols,
+          unsigned int   const rows,
+          unsigned int   const n,
+          unsigned int   const rseed,
+          struct Gauss * const gaussP) {
 /*----------------------------------------------------------------------------
   Generate planet from elevation array.
 
@@ -971,13 +1006,13 @@ genplanet(bool         const stars,
                    clouds ? "clouds" : "planet",
                    rseed, fracdim, powscale, meshsize);
         createPlanetStuff(clouds, a, n, &u, &u1, &bxf, &bxc, &cp, &sunvec,
-                          cols, maxval);
+                          cols, maxval, gaussP);
     }
 
     pixelrow = ppm_allocrow(cols);
     for (row = 0; row < rows; ++row) {
         if (stars)
-            generateStarrySkyRow(pixelrow, cols);
+            generateStarrySkyRow(pixelrow, cols, gaussP);
         else {
             double const by = (n - 1) * uprj(row, rows);
             int    const byf = floor(by) * n;
@@ -993,7 +1028,8 @@ genplanet(bool         const stars,
                 generatePlanetRow(pixelrow, row, rows, cols,
                                   t, t1, u, u1, cp, bxc, bxf, byc, byf,
                                   sunvec,
-                                  maxval);
+                                  maxval,
+                                  gaussP);
         }
         ppm_writeppmrow(stdout, pixelrow, cols, maxval, FALSE);
     }
@@ -1098,14 +1134,15 @@ planet(unsigned int const cols,
 -----------------------------------------------------------------------------*/
     float * a;
     bool error;
+    struct Gauss gauss;
 
-    initgauss(rseed);
+    initgauss(&gauss, rseed);
 
     if (stars) {
         a = NULL;
         error = FALSE;
     } else {
-        spectralsynth(&a, meshsize, 3.0 - fracdim);
+        spectralsynth(&a, meshsize, 3.0 - fracdim, &gauss);
         if (a == NULL) {
             error = TRUE;
         } else {
@@ -1117,7 +1154,7 @@ planet(unsigned int const cols,
         }
     }
     if (!error)
-        genplanet(stars, clouds, a, cols, rows, meshsize, rseed);
+        genplanet(stars, clouds, a, cols, rows, meshsize, rseed, &gauss);
 
     if (a != NULL)
         free(a);
@@ -1159,10 +1196,10 @@ main(int argc, const char ** argv) {
     cols = (MAX(cmdline.height, cmdline.width) + 1) & (~1);
     rows = cmdline.height;
 
-    success = planet(cols, rows, cmdline.night, cmdline.clouds, cmdline.seed);
+    success = planet(cols, rows, cmdline.night,
+                     cmdline.clouds, cmdline.seed);
 
     exit(success ? 0 : 1);
 }
 
 
-
diff --git a/generator/ppmpat.c b/generator/ppmpat.c
index 908c200f..2fd2088d 100644
--- a/generator/ppmpat.c
+++ b/generator/ppmpat.c
@@ -23,8 +23,9 @@
 
 #include "pm_c_util.h"
 #include "mallocvar.h"
-#include "shhopt.h"
 #include "nstring.h"
+#include "rand.h"
+#include "shhopt.h"
 #include "ppm.h"
 #include "ppmdraw.h"
 
@@ -223,9 +224,9 @@ parseCommandLine(int argc, const char ** argv,
     OPTENT3(0, "spiro1",        OPT_FLAG,   NULL,
             &spiro1,     0);
     OPTENT3(0, "spiro2",        OPT_FLAG,   NULL,
-            &spiro1,     0);
+            &spiro2,     0);
     OPTENT3(0, "spiro3",        OPT_FLAG,   NULL,
-            &spiro1,     0);
+            &spiro3,     0);
 #else
     spiro1 = spiro2 = spiro3 = 0;
 #endif
@@ -339,14 +340,15 @@ validateComputableDimensions(unsigned int const cols,
 
 
 static pixel
-randomColor(pixval const maxval) {
+randomColor(struct pm_randSt * const randStP,
+            pixval             const maxval) {
 
     pixel p;
 
     PPM_ASSIGN(p,
-               rand() % (maxval + 1),
-               rand() % (maxval + 1),
-               rand() % (maxval + 1)
+               pm_rand(randStP) % (maxval + 1),
+               pm_rand(randStP) % (maxval + 1),
+               pm_rand(randStP) % (maxval + 1)
         );
 
     return p;
@@ -354,15 +356,18 @@ randomColor(pixval const maxval) {
 
 
 
-#define DARK_THRESH 0.25
+static double const DARK_THRESH = 0.25;
+
+
 
 static pixel
-randomBrightColor(pixval const maxval) {
+randomBrightColor(struct pm_randSt * const randStP,
+                  pixval             const maxval) {
 
     pixel p;
 
     do {
-        p = randomColor(maxval);
+        p = randomColor(randStP, maxval);
     } while (PPM_LUMIN(p) <= maxval * DARK_THRESH);
 
     return p;
@@ -371,12 +376,13 @@ randomBrightColor(pixval const maxval) {
 
 
 static pixel
-randomDarkColor(pixval const maxval) {
+randomDarkColor(struct pm_randSt * const randStP,
+                pixval             const maxval) {
 
     pixel p;
 
     do {
-        p = randomColor(maxval);
+        p = randomColor(randStP, maxval);
     } while (PPM_LUMIN(p) > maxval * DARK_THRESH);
 
     return p;
@@ -448,7 +454,8 @@ nextColorBg(ColorTable * const colorTableP) {
 
 
 static pixel
-randomAnticamoColor(pixval const maxval) {
+randomAnticamoColor(struct pm_randSt * const randStP,
+                    pixval             const maxval) {
 
     int v1, v2, v3;
     pixel p;
@@ -457,37 +464,49 @@ randomAnticamoColor(pixval const maxval) {
     v2 = (maxval + 1) / 2;
     v3 = 3 * v1;
 
-    switch (rand() % 15) {
+    switch (pm_rand(randStP) % 15) {
     case 0: case 1:
-        PPM_ASSIGN(p, rand() % v1 + v3, rand() % v2, rand() % v2);
+        PPM_ASSIGN(p, pm_rand(randStP) % v1 + v3,
+                      pm_rand(randStP) % v2,
+                      pm_rand(randStP) % v2);
         break;
 
     case 2:
     case 3:
-        PPM_ASSIGN(p, rand() % v2, rand() % v1 + v3, rand() % v2);
+        PPM_ASSIGN(p, pm_rand(randStP) % v2,
+                      pm_rand(randStP) % v1 + v3,
+                      pm_rand(randStP) % v2);
         break;
 
     case 4:
     case 5:
-        PPM_ASSIGN(p, rand() % v2, rand() % v2, rand() % v1 + v3);
+        PPM_ASSIGN(p, pm_rand(randStP) % v2,
+                      pm_rand(randStP) % v2,
+                      pm_rand(randStP) % v1 + v3);
         break;
 
     case 6:
     case 7:
     case 8:
-        PPM_ASSIGN(p, rand() % v2, rand() % v1 + v3, rand() % v1 + v3);
+        PPM_ASSIGN(p, pm_rand(randStP) % v2,
+                      pm_rand(randStP) % v1 + v3,
+                      pm_rand(randStP) % v1 + v3);
         break;
 
     case 9:
     case 10:
     case 11:
-        PPM_ASSIGN(p, rand() % v1 + v3, rand() % v2, rand() % v1 + v3);
+        PPM_ASSIGN(p, pm_rand(randStP) % v1 + v3,
+                      pm_rand(randStP) % v2,
+                      pm_rand(randStP) % v1 + v3);
         break;
 
     case 12:
     case 13:
     case 14:
-        PPM_ASSIGN(p, rand() % v1 + v3, rand() % v1 + v3, rand() % v2);
+        PPM_ASSIGN(p, pm_rand(randStP) % v1 + v3,
+                      pm_rand(randStP) % v1 + v3,
+                      pm_rand(randStP) % v2);
         break;
     }
 
@@ -497,7 +516,8 @@ randomAnticamoColor(pixval const maxval) {
 
 
 static pixel
-randomCamoColor(pixval const maxval) {
+randomCamoColor(struct pm_randSt * const randStP,
+                pixval             const maxval) {
 
     int const v1 = (maxval + 1 ) / 8;
     int const v2 = (maxval + 1 ) / 4;
@@ -505,31 +525,39 @@ randomCamoColor(pixval const maxval) {
 
     pixel p;
 
-    switch (rand() % 10) {
+    switch (pm_rand(randStP) % 10) {
     case 0:
     case 1:
     case 2:
         /* light brown */
-        PPM_ASSIGN(p, rand() % v3 + v3, rand() % v3 + v2, rand() % v3 + v2);
+        PPM_ASSIGN(p, pm_rand(randStP) % v3 + v3,
+                      pm_rand(randStP) % v3 + v2,
+                      pm_rand(randStP) % v3 + v2);
         break;
 
     case 3:
     case 4:
     case 5:
         /* dark green */
-        PPM_ASSIGN(p, rand() % v2, rand() % v2 + 3 * v1, rand() % v2);
+        PPM_ASSIGN(p, pm_rand(randStP) % v2,
+                      pm_rand(randStP) % v2 + 3 * v1,
+                      pm_rand(randStP) % v2);
         break;
 
     case 6:
     case 7:
         /* brown */
-        PPM_ASSIGN(p, rand() % v2 + v2, rand() % v2, rand() % v2);
+        PPM_ASSIGN(p, pm_rand(randStP) % v2 + v2,
+                      pm_rand(randStP) % v2,
+                      pm_rand(randStP) % v2);
         break;
 
     case 8:
     case 9:
         /* dark brown */
-        PPM_ASSIGN(p, rand() % v1 + v1, rand() % v1, rand() % v1);
+        PPM_ASSIGN(p, pm_rand(randStP) % v1 + v1,
+                      pm_rand(randStP) % v1,
+                      pm_rand(randStP) % v1);
         break;
     }
 
@@ -539,28 +567,30 @@ randomCamoColor(pixval const maxval) {
 
 
 static float
-rnduni(void) {
-    return rand() % 32767 / 32767.0;
+rnduni(struct pm_randSt * const randStP) {
+
+    return pm_rand(randStP) % 32767 / 32767.0;
 }
 
 
 
 static void
-clearBackgroundCamo(pixel **     const pixels,
-                    unsigned int const cols,
-                    unsigned int const rows,
-                    pixval       const maxval,
-                    ColorTable * const colorTableP,
-                    bool         const antiflag) {
+clearBackgroundCamo(pixel **           const pixels,
+                    unsigned int       const cols,
+                    unsigned int       const rows,
+                    pixval             const maxval,
+                    ColorTable *       const colorTableP,
+                    struct pm_randSt * const randStP,
+                    bool               const antiflag) {
 
     pixel color;
 
     if (colorTableP->count > 0) {
         color = colorTableP->color[0];
     } else if (antiflag)
-        color = randomAnticamoColor(maxval);
+        color = randomAnticamoColor(randStP, maxval);
     else
-        color = randomCamoColor(maxval);
+        color = randomCamoColor(randStP,maxval);
 
     ppmd_filledrectangle(
         pixels, cols, rows, maxval, 0, 0, cols, rows, PPMD_NULLDRAWPROC,
@@ -570,13 +600,14 @@ clearBackgroundCamo(pixel **     const pixels,
 
 
 static void
-camoFill(pixel **         const pixels,
-         unsigned int     const cols,
-         unsigned int     const rows,
-         pixval           const maxval,
-         struct fillobj * const fh,
-         ColorTable     * const colorTableP,
-         bool             const antiflag) {
+camoFill(pixel **           const pixels,
+         unsigned int       const cols,
+         unsigned int       const rows,
+         pixval             const maxval,
+         struct fillobj *   const fh,
+         ColorTable     *   const colorTableP,
+         struct pm_randSt * const randStP,
+         bool               const antiflag) {
 
     pixel color;
 
@@ -585,9 +616,9 @@ camoFill(pixel **         const pixels,
         color = colorTableP->color[colorTableP->index];
         nextColorBg(colorTableP);
     } else if (antiflag)
-        color = randomAnticamoColor(maxval);
+        color = randomAnticamoColor(randStP, maxval);
     else
-        color = randomCamoColor(maxval);
+        color = randomCamoColor(randStP, maxval);
 
     ppmd_fill(pixels, cols, rows, maxval, fh, PPMD_NULLDRAWPROC, &color);
 }
@@ -608,25 +639,29 @@ camoFill(pixel **         const pixels,
 
 
 static void
-computeXsYs(int *        const xs,
-            int *        const ys,
-            unsigned int const cols,
-            unsigned int const rows,
-            unsigned int const pointCt) {
-
-    unsigned int const cx = rand() % cols;
-    unsigned int const cy = rand() % rows;
-    double const a = rnduni() * (MAX_ELLIPSE_FACTOR - MIN_ELLIPSE_FACTOR) +
-        MIN_ELLIPSE_FACTOR;
-    double const b = rnduni() * (MAX_ELLIPSE_FACTOR - MIN_ELLIPSE_FACTOR) +
-        MIN_ELLIPSE_FACTOR;
-    double const theta = rnduni() * 2.0 * M_PI;
+computeXsYs(int *              const xs,
+            int *              const ys,
+            unsigned int       const cols,
+            unsigned int       const rows,
+            unsigned int       const pointCt,
+            struct pm_randSt * const randStP) {
+
+    unsigned int const cx = pm_rand(randStP) % cols;
+    unsigned int const cy = pm_rand(randStP) % rows;
+    double const a = rnduni(randStP) *
+                       (MAX_ELLIPSE_FACTOR - MIN_ELLIPSE_FACTOR) +
+                        MIN_ELLIPSE_FACTOR;
+    double const b = rnduni(randStP) *
+                       (MAX_ELLIPSE_FACTOR - MIN_ELLIPSE_FACTOR) +
+                        MIN_ELLIPSE_FACTOR;
+    double const theta = rnduni(randStP) * 2.0 * M_PI;
 
     unsigned int p;
 
     for (p = 0; p < pointCt; ++p) {
-        double const c = rnduni() * (MAX_POINT_FACTOR - MIN_POINT_FACTOR) +
-            MIN_POINT_FACTOR;
+        double const c = rnduni(randStP) *
+                           (MAX_POINT_FACTOR - MIN_POINT_FACTOR) +
+                            MIN_POINT_FACTOR;
         double const tx   = a * sin(p * 2.0 * M_PI / pointCt);
         double const ty   = b * cos(p * 2.0 * M_PI / pointCt);
         double const tang = atan2(ty, tx) + theta;
@@ -638,18 +673,20 @@ computeXsYs(int *        const xs,
 
 
 static void
-camo(pixel **     const pixels,
-     unsigned int const cols,
-     unsigned int const rows,
-     ColorTable * const colorTableP,
-     pixval       const maxval,
-     bool         const antiflag) {
+camo(pixel **           const pixels,
+     unsigned int       const cols,
+     unsigned int       const rows,
+     ColorTable *       const colorTableP,
+     struct pm_randSt * const randStP,
+     pixval             const maxval,
+     bool               const antiflag) {
 
     unsigned int const n = (rows * cols) / SQR(BLOBRAD) * 5;
 
     unsigned int i;
 
-    clearBackgroundCamo(pixels, cols, rows, maxval, colorTableP, antiflag);
+    clearBackgroundCamo(pixels, cols, rows, maxval,
+                        colorTableP, randStP, antiflag);
 
     if (colorTableP->count > 0) {
         assert(colorTableP->count > 1);
@@ -658,13 +695,13 @@ camo(pixel **     const pixels,
 
     for (i = 0; i < n; ++i) {
         unsigned int const pointCt =
-            rand() % (MAX_POINTS - MIN_POINTS + 1) + MIN_POINTS;
+            pm_rand(randStP) % (MAX_POINTS - MIN_POINTS + 1) + MIN_POINTS;
 
         int xs[MAX_POINTS], ys[MAX_POINTS];
         int x0, y0;
         struct fillobj * fh;
 
-        computeXsYs(xs, ys, cols, rows, pointCt);
+        computeXsYs(xs, ys, cols, rows, pointCt, randStP);
 
         x0 = (xs[0] + xs[pointCt - 1]) / 2;
         y0 = (ys[0] + ys[pointCt - 1]) / 2;
@@ -675,7 +712,8 @@ camo(pixel **     const pixels,
             pixels, cols, rows, maxval, x0, y0, pointCt, xs, ys, x0, y0,
             ppmd_fill_drawproc, fh);
 
-        camoFill(pixels, cols, rows, maxval, fh, colorTableP, antiflag);
+        camoFill(pixels, cols, rows, maxval, fh,
+                 colorTableP, randStP, antiflag);
 
         ppmd_fill_destroy(fh);
     }
@@ -688,17 +726,20 @@ camo(pixel **     const pixels,
 -----------------------------------------------------------------------------*/
 
 static void
-gingham2(pixel **     const pixels,
-         unsigned int const cols,
-         unsigned int const rows,
-         ColorTable   const colorTable,
-         pixval       const maxval) {
+gingham2(pixel **           const pixels,
+         unsigned int       const cols,
+         unsigned int       const rows,
+         ColorTable         const colorTable,
+         struct pm_randSt * const randStP,
+         pixval             const maxval) {
 
     bool  const colorSpec = (colorTable.count > 0);
     pixel const backcolor = colorSpec ?
-                            colorTable.color[0] : randomDarkColor(maxval);
+                              colorTable.color[0] :
+                              randomDarkColor(randStP, maxval);
     pixel const forecolor = colorSpec ?
-                            colorTable.color[1] : randomBrightColor(maxval);
+                              colorTable.color[1] :
+                              randomBrightColor(randStP, maxval);
     unsigned int const colso2 = cols / 2;
     unsigned int const rowso2 = rows / 2;
 
@@ -722,19 +763,23 @@ gingham2(pixel **     const pixels,
 
 
 static void
-gingham3(pixel **     const pixels,
-         unsigned int const cols,
-         unsigned int const rows,
-         ColorTable   const colorTable,
-         pixval       const maxval) {
+gingham3(pixel **           const pixels,
+         unsigned int       const cols,
+         unsigned int       const rows,
+         ColorTable         const colorTable,
+         struct pm_randSt * const randStP,
+         pixval             const maxval) {
 
     bool  const colorSpec = (colorTable.count > 0);
-    pixel const backcolor = colorSpec ?
-                            colorTable.color[0] : randomDarkColor(maxval);
+    pixel const backcolor  = colorSpec ?
+                               colorTable.color[0] :
+                               randomDarkColor(randStP, maxval);
     pixel const fore1color = colorSpec ?
-                            colorTable.color[1] : randomBrightColor(maxval);
+                               colorTable.color[1] :
+                               randomBrightColor(randStP, maxval);
     pixel const fore2color = colorSpec ?
-                            colorTable.color[2] : randomBrightColor(maxval);
+                               colorTable.color[2] :
+                               randomBrightColor(randStP, maxval);
     unsigned int const colso4 = cols / 4;
     unsigned int const rowso4 = rows / 4;
 
@@ -770,19 +815,23 @@ gingham3(pixel **     const pixels,
 
 
 static void
-madras(pixel **     const pixels,
-       unsigned int const cols,
-       unsigned int const rows,
-       ColorTable   const colorTable,
-       pixval       const maxval) {
+madras(pixel **           const pixels,
+       unsigned int       const cols,
+       unsigned int       const rows,
+       ColorTable         const colorTable,
+       struct pm_randSt * const randStP,
+       pixval             const maxval) {
 
     bool  const colorSpec = (colorTable.count > 0);
-    pixel const backcolor = colorSpec ?
-                            colorTable.color[0] : randomDarkColor(maxval);
+    pixel const backcolor  = colorSpec ?
+                               colorTable.color[0] :
+                               randomDarkColor(randStP, maxval);
     pixel const fore1color = colorSpec ?
-                            colorTable.color[1] : randomBrightColor(maxval);
+                               colorTable.color[1] :
+                               randomBrightColor(randStP, maxval);
     pixel const fore2color = colorSpec ?
-                            colorTable.color[2] : randomBrightColor(maxval);
+                               colorTable.color[2] :
+                               randomBrightColor(randStP, maxval);
     unsigned int const cols2  = cols * 2 / 44;
     unsigned int const rows2  = rows * 2 / 44;
     unsigned int const cols3  = cols * 3 / 44;
@@ -899,19 +948,23 @@ madras(pixel **     const pixels,
 
 
 static void
-tartan(pixel **     const pixels,
-       unsigned int const cols,
-       unsigned int const rows,
-       ColorTable   const colorTable,
-       pixval       const maxval) {
+tartan(pixel **           const pixels,
+       unsigned int       const cols,
+       unsigned int       const rows,
+       ColorTable         const colorTable,
+       struct pm_randSt * const randStP,
+       pixval             const maxval) {
 
     bool  const colorSpec = (colorTable.count > 0);
-    pixel const backcolor = colorSpec ?
-                            colorTable.color[0] : randomDarkColor(maxval);
+    pixel const backcolor  = colorSpec ?
+                               colorTable.color[0] :
+                               randomDarkColor(randStP, maxval);
     pixel const fore1color = colorSpec ?
-                            colorTable.color[1] : randomBrightColor(maxval);
+                               colorTable.color[1] :
+                               randomBrightColor(randStP, maxval);
     pixel const fore2color = colorSpec ?
-                            colorTable.color[2] : randomBrightColor(maxval);
+                               colorTable.color[2] :
+                               randomBrightColor(randStP, maxval);
     unsigned int const cols1  = cols / 22;
     unsigned int const rows1  = rows / 22;
     unsigned int const cols3  = cols * 3 / 22;
@@ -1013,14 +1066,15 @@ argyle(pixel **     const pixels,
        unsigned int const cols,
        unsigned int const rows,
        ColorTable   const colorTable,
+       struct pm_randSt * const randStP,
        pixval       const maxval,
        bool         const stripes) {
 
     bool  const colorSpec = (colorTable.count > 0);
     pixel const backcolor = colorSpec ?
-        colorTable.color[0] : randomDarkColor(maxval);
+        colorTable.color[0] : randomDarkColor(randStP, maxval);
     pixel const forecolor = colorSpec ?
-        colorTable.color[1] : randomBrightColor(maxval);
+        colorTable.color[1] : randomBrightColor(randStP, maxval);
 
     /* Fill canvas with background to start */
     ppmd_filledrectangle(
@@ -1032,7 +1086,8 @@ argyle(pixel **     const pixels,
     if (stripes) {
          /* Connect corners with thin stripes */
          pixel const stripecolor =
-             colorSpec ? colorTable.color[2] : randomBrightColor(maxval);
+             colorSpec ? colorTable.color[2] :
+                         randomBrightColor(randStP, maxval);
 
          ppmd_line(pixels, cols, rows, maxval, 0, 0, cols-1, rows-1,
               PPMD_NULLDRAWPROC, (char *) &stripecolor);
@@ -1049,32 +1104,33 @@ argyle(pixel **     const pixels,
 
 
 
-#define MAXPOLES 500
+static unsigned int const MAXPOLES = 500;
 
 
 
 static void
-placeAndColorPolesRandomly(int *        const xs,
-                           int *        const ys,
-                           pixel *      const colors,
-                           unsigned int const cols,
-                           unsigned int const rows,
-                           pixval       const maxval,
-                           ColorTable * const colorTableP,
-                           unsigned int const poleCt) {
+placeAndColorPolesRandomly(int *              const xs,
+                           int *              const ys,
+                           pixel *            const colors,
+                           unsigned int       const cols,
+                           unsigned int       const rows,
+                           pixval             const maxval,
+                           ColorTable *       const colorTableP,
+                           struct pm_randSt * const randStP,
+                           unsigned int       const poleCt) {
 
     unsigned int i;
 
     for (i = 0; i < poleCt; ++i) {
 
-        xs[i] = rand() % cols;
-        ys[i] = rand() % rows;
+        xs[i] = pm_rand(randStP) % cols;
+        ys[i] = pm_rand(randStP) % rows;
 
         if (colorTableP->count > 0) {
             colors[i] = colorTableP->color[colorTableP->index];
             nextColor(colorTableP);
         } else
-            colors[i] = randomBrightColor(maxval);
+            colors[i] = randomBrightColor(randStP, maxval);
     }
 }
 
@@ -1104,11 +1160,12 @@ assignInterpolatedColor(pixel * const resultP,
 
 
 static void
-poles(pixel **     const pixels,
-      unsigned int const cols,
-      unsigned int const rows,
-      ColorTable * const colorTableP,
-      pixval       const maxval) {
+poles(pixel **           const pixels,
+      unsigned int       const cols,
+      unsigned int       const rows,
+      ColorTable *       const colorTableP,
+      struct pm_randSt * const randStP,
+      pixval             const maxval) {
 
     unsigned int const poleCt = MAX(2, MIN(MAXPOLES, cols * rows / 30000));
 
@@ -1117,7 +1174,7 @@ poles(pixel **     const pixels,
     unsigned int row;
 
     placeAndColorPolesRandomly(xs, ys, colors, cols, rows, maxval,
-                               colorTableP, poleCt);
+                               colorTableP, randStP, poleCt);
 
     /* Interpolate points */
 
@@ -1236,11 +1293,12 @@ sqRainbowCircleDrawproc(pixel **     const pixels,
 
 
 static void
-chooseSqPoleColors(ColorTable * const colorTableP,
-                   pixval       const maxval,
-                   pixel *      const color1P,
-                   pixel *      const color2P,
-                   pixel *      const color3P) {
+chooseSqPoleColors(ColorTable *       const colorTableP,
+                   pixval             const maxval,
+                   pixel *            const color1P,
+                   pixel *            const color2P,
+                   pixel *            const color3P,
+                   struct pm_randSt * const randStP) {
 
     if (colorTableP->count > 0) {
         *color1P = colorTableP->color[colorTableP->index];
@@ -1250,19 +1308,20 @@ chooseSqPoleColors(ColorTable * const colorTableP,
         *color3P = colorTableP->color[colorTableP->index];
         nextColor(colorTableP);
     } else {
-        *color1P = randomBrightColor(maxval);
-        *color2P = randomBrightColor(maxval);
-        *color3P = randomBrightColor(maxval);
+        *color1P = randomBrightColor(randStP, maxval);
+        *color2P = randomBrightColor(randStP, maxval);
+        *color3P = randomBrightColor(randStP, maxval);
     }
 }
 
 
 
 static void
-sqAssignColors(unsigned int const circlecount,
-               pixval       const maxval,
-               ColorTable * const colorTableP,
-               pixel *      const colors) {
+sqAssignColors(unsigned int       const circlecount,
+               pixval             const maxval,
+               ColorTable *       const colorTableP,
+               pixel *            const colors,
+               struct pm_randSt * const randStP) {
 
     float const cco3 = (circlecount - 1) / 3.0;
 
@@ -1271,7 +1330,7 @@ sqAssignColors(unsigned int const circlecount,
     pixel rc3;
     unsigned int i;
 
-    chooseSqPoleColors(colorTableP, maxval, &rc1, &rc2, &rc3);
+    chooseSqPoleColors(colorTableP, maxval, &rc1, &rc2, &rc3, randStP);
 
     for (i = 0; i < circlecount; ++i) {
         if (i < cco3) {
@@ -1333,24 +1392,25 @@ clearBackgroundSquig(pixel **     const pixels,
 
 
 static void
-chooseWrapAroundPoint(unsigned int const cols,
-                      unsigned int const rows,
-                      ppmd_point * const pFirstP,
-                      ppmd_point * const pLastP,
-                      ppmd_point * const p0P,
-                      ppmd_point * const p1P,
-                      ppmd_point * const p2P,
-                      ppmd_point * const p3P) {
-
-    switch (rand() % 4) {
+chooseWrapAroundPoint(unsigned int       const cols,
+                      unsigned int       const rows,
+                      ppmd_point *       const pFirstP,
+                      ppmd_point *       const pLastP,
+                      ppmd_point *       const p0P,
+                      ppmd_point *       const p1P,
+                      ppmd_point *       const p2P,
+                      ppmd_point *       const p3P,
+                      struct pm_randSt * const randStP) {
+
+    switch (pm_rand(randStP) % 4) {
     case 0:
-        p1P->x = rand() % cols;
+        p1P->x = pm_rand(randStP) % cols;
         p1P->y = 0;
         if (p1P->x < cols / 2)
-            pFirstP->x = rand() % (p1P->x * 2 + 1);
+            pFirstP->x = pm_rand(randStP) % (p1P->x * 2 + 1);
         else
-            pFirstP->x = cols - 1 - rand() % ((cols - p1P->x) * 2);
-        pFirstP->y = rand() % rows;
+            pFirstP->x = cols - 1 - pm_rand(randStP) % ((cols - p1P->x) * 2);
+        pFirstP->y = pm_rand(randStP) % rows;
         p2P->x = p1P->x;
         p2P->y = rows - 1;
         pLastP->x = 2 * p2P->x - pFirstP->x;
@@ -1362,13 +1422,13 @@ chooseWrapAroundPoint(unsigned int const cols,
         break;
 
     case 1:
-        p2P->x = rand() % cols;
+        p2P->x = pm_rand(randStP) % cols;
         p2P->y = 0;
         if (p2P->x < cols / 2)
-            pLastP->x = rand() % (p2P->x * 2 + 1);
+            pLastP->x = pm_rand(randStP) % (p2P->x * 2 + 1);
         else
-            pLastP->x = cols - 1 - rand() % ((cols - p2P->x) * 2);
-        pLastP->y = rand() % rows;
+            pLastP->x = cols - 1 - pm_rand(randStP) % ((cols - p2P->x) * 2);
+        pLastP->y = pm_rand(randStP) % rows;
         p1P->x = p2P->x;
         p1P->y = rows - 1;
         pFirstP->x = 2 * p1P->x - pLastP->x;
@@ -1381,12 +1441,12 @@ chooseWrapAroundPoint(unsigned int const cols,
 
     case 2:
         p1P->x = 0;
-        p1P->y = rand() % rows;
-        pFirstP->x = rand() % cols;
+        p1P->y = pm_rand(randStP) % rows;
+        pFirstP->x = pm_rand(randStP) % cols;
         if (p1P->y < rows / 2)
-            pFirstP->y = rand() % (p1P->y * 2 + 1);
+            pFirstP->y = pm_rand(randStP) % (p1P->y * 2 + 1);
         else
-            pFirstP->y = rows - 1 - rand() % ((rows - p1P->y) * 2);
+            pFirstP->y = rows - 1 - pm_rand(randStP) % ((rows - p1P->y) * 2);
         p2P->x = cols - 1;
         p2P->y = p1P->y;
         pLastP->x = p2P->x - pFirstP->x;
@@ -1399,12 +1459,12 @@ chooseWrapAroundPoint(unsigned int const cols,
 
     case 3:
         p2P->x = 0;
-        p2P->y = rand() % rows;
-        pLastP->x = rand() % cols;
+        p2P->y = pm_rand(randStP) % rows;
+        pLastP->x = pm_rand(randStP) % cols;
         if (p2P->y < rows / 2)
-            pLastP->y = rand() % (p2P->y * 2 + 1);
+            pLastP->y = pm_rand(randStP) % (p2P->y * 2 + 1);
         else
-            pLastP->y = rows - 1 - rand() % ((rows - p2P->y) * 2);
+            pLastP->y = rows - 1 - pm_rand(randStP) % ((rows - p2P->y) * 2);
         p1P->x = cols - 1;
         p1P->y = p2P->y;
         pFirstP->x = p1P->x - pLastP->x;
@@ -1420,11 +1480,12 @@ chooseWrapAroundPoint(unsigned int const cols,
 
 
 static void
-squig(pixel **     const pixels,
-      unsigned int const cols,
-      unsigned int const rows,
-      ColorTable * const colorTableP,
-      pixval       const maxval) {
+squig(pixel **           const pixels,
+      unsigned int       const cols,
+      unsigned int       const rows,
+      ColorTable *       const colorTableP,
+      struct pm_randSt * const randStP,
+      pixval             const maxval) {
 
     int i;
 
@@ -1453,10 +1514,11 @@ squig(pixel **     const pixels,
         ppmd_circlep(pixels, cols, rows, maxval,
                      ppmd_makePoint(0, 0), radius,
                      sqMeasureCircleDrawproc, &sqClientData);
-        sqAssignColors(squig.circleCt, maxval, colorTableP, squig.color);
+        sqAssignColors(squig.circleCt, maxval, colorTableP, squig.color,
+                       randStP);
 
         chooseWrapAroundPoint(cols, rows, &c[0], &c[SQ_POINTS-1],
-                              &p0, &p1, &p2, &p3);
+                              &p0, &p1, &p2, &p3, randStP);
 
         {
             /* Do the middle points */
@@ -1466,8 +1528,8 @@ squig(pixel **     const pixels,
               /* validateSquigAspect() assures that
                  cols - 2 * radius, rows -2 * radius are positive
               */
-                c[j].x = (rand() % (cols - 2 * radius)) + radius;
-                c[j].y = (rand() % (rows - 2 * radius)) + radius;
+                c[j].x = (pm_rand(randStP) % (cols - 2 * radius)) + radius;
+                c[j].y = (pm_rand(randStP) % (rows - 2 * radius)) + radius;
             }
         }
 
@@ -1490,6 +1552,7 @@ main(int argc, const char ** argv) {
 
     struct CmdlineInfo cmdline;
     pixel ** pixels;
+    struct pm_randSt randSt;
 
     pm_proginit(&argc, argv);
 
@@ -1497,65 +1560,68 @@ main(int argc, const char ** argv) {
 
     validateComputableDimensions(cmdline.width, cmdline.height);
 
-    srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
+    pm_randinit(&randSt);
+    pm_srand2(&randSt, cmdline.randomseedSpec, cmdline.randomseed);
 
     pixels = ppm_allocarray(cmdline.width, cmdline.height);
 
     switch (cmdline.basePattern) {
     case PAT_GINGHAM2:
         gingham2(pixels, cmdline.width, cmdline.height,
-                 cmdline.colorTable, PPM_MAXMAXVAL);
+                 cmdline.colorTable, &randSt, PPM_MAXMAXVAL);
         break;
 
     case PAT_GINGHAM3:
         gingham3(pixels, cmdline.width, cmdline.height,
-                 cmdline.colorTable, PPM_MAXMAXVAL);
+                 cmdline.colorTable, &randSt, PPM_MAXMAXVAL);
         break;
 
     case PAT_MADRAS:
         madras(pixels, cmdline.width, cmdline.height,
-               cmdline.colorTable, PPM_MAXMAXVAL);
+               cmdline.colorTable, &randSt, PPM_MAXMAXVAL);
         break;
 
     case PAT_TARTAN:
         tartan(pixels, cmdline.width, cmdline.height,
-               cmdline.colorTable, PPM_MAXMAXVAL);
+               cmdline.colorTable, &randSt, PPM_MAXMAXVAL);
         break;
 
     case PAT_ARGYLE1:
         argyle(pixels, cmdline.width, cmdline.height,
-               cmdline.colorTable, PPM_MAXMAXVAL, FALSE);
+               cmdline.colorTable, &randSt, PPM_MAXMAXVAL, FALSE);
         break;
 
     case PAT_ARGYLE2:
         argyle(pixels, cmdline.width, cmdline.height,
-               cmdline.colorTable, PPM_MAXMAXVAL, TRUE);
+               cmdline.colorTable, &randSt, PPM_MAXMAXVAL, TRUE);
         break;
 
     case PAT_POLES:
         poles(pixels, cmdline.width, cmdline.height,
-              &cmdline.colorTable, PPM_MAXMAXVAL);
+              &cmdline.colorTable, &randSt, PPM_MAXMAXVAL);
         break;
 
     case PAT_SQUIG:
         squig(pixels, cmdline.width, cmdline.height,
-              &cmdline.colorTable, PPM_MAXMAXVAL);
+              &cmdline.colorTable, &randSt, PPM_MAXMAXVAL);
         break;
 
     case PAT_CAMO:
         camo(pixels, cmdline.width, cmdline.height,
-             &cmdline.colorTable, PPM_MAXMAXVAL, 0);
+             &cmdline.colorTable, &randSt, PPM_MAXMAXVAL, 0);
         break;
 
     case PAT_ANTICAMO:
         camo(pixels, cmdline.width, cmdline.height,
-             &cmdline.colorTable, PPM_MAXMAXVAL, 1);
+             &cmdline.colorTable, &randSt, PPM_MAXMAXVAL, 1);
         break;
 
     default:
         pm_error("can't happen!");
     }
 
+    pm_randterm(&randSt);
+
     ppm_writeppm(stdout, pixels, cmdline.width, cmdline.height,
                  PPM_MAXMAXVAL, 0);
 
@@ -1567,4 +1633,3 @@ main(int argc, const char ** argv) {
 }
 
 
-
diff --git a/generator/ppmrough.c b/generator/ppmrough.c
index c87a0364..88a27a69 100644
--- a/generator/ppmrough.c
+++ b/generator/ppmrough.c
@@ -16,10 +16,10 @@
 
 #include "pm_c_util.h"
 #include "mallocvar.h"
+#include "rand.h"
 #include "shhopt.h"
 #include "ppm.h"
 
-static pixel** PIX;
 static pixval BG_RED, BG_GREEN, BG_BLUE;
 
 
@@ -89,95 +89,239 @@ parseCommandLine(int argc, const char ** argv,
 
 
 static void
-procLeft(int          const r1,
-         int          const r2,
-         int          const c1,
-         int          const c2,
-         unsigned int const var) {
+makeAllForegroundColor(pixel **     const pixels,
+                       unsigned int const rows,
+                       unsigned int const cols,
+                       pixval       const r,
+                       pixval       const g,
+                       pixval       const b) {
 
-    int cm, rm, c;
+    unsigned int row;
 
-    if (r1 + 1 == r2)  return;
-    rm = (r1 + r2) >> 1;
-    cm = (c1 + c2) >> 1;
-    cm += (int)floor(((float)rand() / RAND_MAX - 0.5) * var + 0.5);
+    for (row = 0; row < rows; ++row) {
+        unsigned int col;
+
+        for (col = 0; col < cols; ++col)
+            PPM_ASSIGN(pixels[row][col], r, g, b);
+    }
+}
 
-    for (c = 0; c < cm; c++)
-        PPM_ASSIGN(PIX[rm][c], BG_RED, BG_GREEN, BG_BLUE);
 
-    procLeft(r1, rm, c1, cm, var);
-    procLeft(rm, r2, cm, c2, var);
+
+static void
+procLeft(pixel **           const pixels,
+         int                const r1,
+         int                const r2,
+         int                const c1,
+         int                const c2,
+         unsigned int       const var,
+         struct pm_randSt * const randStP) {
+
+    if (r1 + 1 != r2) {
+        int const rm = (r1 + r2) >> 1;
+        int const cm = ((c1 + c2) >> 1) +
+            (int)floor(((float)pm_rand(randStP) / RAND_MAX - 0.5) * var + 0.5);
+
+        int c;
+
+        for (c = 0; c < cm; c++)
+            PPM_ASSIGN(pixels[rm][c], BG_RED, BG_GREEN, BG_BLUE);
+
+        procLeft(pixels, r1, rm, c1, cm, var, randStP);
+        procLeft(pixels, rm, r2, cm, c2, var, randStP);
+    }
 }
 
 
 
 static void
-procRight(int          const r1,
-          int          const r2,
-          int          const c1,
-          int          const c2,
-          unsigned int const width,
-          unsigned int const var) {
+procRight(pixel **           const pixels,
+          int                const r1,
+          int                const r2,
+          int                const c1,
+          int                const c2,
+          unsigned int       const width,
+          unsigned int       const var,
+          struct pm_randSt * const randStP) {
+
+    if (r1 + 1 != r2) {
+        int const rm = (r1 + r2) >> 1;
+        int const cm = ((c1 + c2) >> 1) +
+            (int)floor(((float)pm_rand(randStP) / RAND_MAX - 0.5) * var + 0.5);
+
+        int c;
+
+        for (c = cm; c < width; c++)
+            PPM_ASSIGN(pixels[rm][c], BG_RED, BG_GREEN, BG_BLUE);
+
+        procRight(pixels, r1, rm, c1, cm, width, var, randStP);
+        procRight(pixels, rm, r2, cm, c2, width, var, randStP);
+    }
+}
 
-    int cm, rm, c;
 
-    if (r1 + 1 == r2)  return;
-    rm = (r1 + r2) >> 1;
-    cm = (c1 + c2) >> 1;
-    cm += (int)floor(((float)rand() / RAND_MAX - 0.5) * var + 0.5);
 
-    for (c = cm; c < width; c++)
-        PPM_ASSIGN(PIX[rm][c], BG_RED, BG_GREEN, BG_BLUE);
+static void
+procTop(pixel **           const pixels,
+        int                const c1,
+        int                const c2,
+        int                const r1,
+        int                const r2,
+        unsigned int       const var,
+        struct pm_randSt * const randStP) {
+
+    if (c1 + 1 != c2) {
+        int const cm = (c1 + c2) >> 1;
+        int const rm = ((r1 + r2) >> 1) +
+            (int)floor(((float)pm_rand(randStP) / RAND_MAX - 0.5) * var + 0.5);
+
+        int r;
+
+        for (r = 0; r < rm; r++)
+            PPM_ASSIGN(pixels[r][cm], BG_RED, BG_GREEN, BG_BLUE);
+
+        procTop(pixels, c1, cm, r1, rm, var, randStP);
+        procTop(pixels, cm, c2, rm, r2, var, randStP);
+    }
+}
+
+
 
-    procRight(r1, rm, c1, cm, width, var);
-    procRight(rm, r2, cm, c2, width, var);
+static void
+procBottom(pixel **           const pixels,
+           int                const c1,
+           int                const c2,
+           int                const r1,
+           int                const r2,
+           unsigned int       const height,
+           unsigned int       const var,
+           struct pm_randSt * const randStP) {
+
+    if (c1 + 1 != c2) {
+        int const cm = (c1 + c2) >> 1;
+        int const rm = ((r1 + r2) >> 1) +
+            (int)floor(((float)pm_rand(randStP) / RAND_MAX - 0.5) * var + 0.5);
+
+        int r;
+
+        for (r = rm; r < height; ++r)
+            PPM_ASSIGN(pixels[r][cm], BG_RED, BG_GREEN, BG_BLUE);
+
+        procBottom(pixels, c1, cm, r1, rm, height, var, randStP);
+        procBottom(pixels, cm, c2, rm, r2, height, var, randStP);
+    }
 }
 
 
 
 static void
-procTop(int          const c1,
-        int          const c2,
-        int          const r1,
-        int          const r2,
-        unsigned int const var) {
+makeRaggedLeftBorder(pixel **           const pixels,
+                     unsigned int       const rows,
+                     unsigned int       const cols,
+                     unsigned int       const left,
+                     unsigned int       const var,
+                     struct pm_randSt * const randStP) {
 
-    int rm, cm, r;
+    if (left >= 0) {
+        int const leftC1 = left;
+        int const leftC2 = left;
+        int const leftR1 = 0;
+        int const leftR2 = rows - 1;
 
-    if (c1 + 1 == c2)  return;
-    cm = (c1 + c2) >> 1;
-    rm = (r1 + r2) >> 1;
-    rm += (int)floor(((float)rand() / RAND_MAX - 0.5) * var + 0.5);
+        unsigned int col;
 
-    for (r = 0; r < rm; r++)
-        PPM_ASSIGN(PIX[r][cm], BG_RED, BG_GREEN, BG_BLUE);
+        for (col = 0; col < leftC1; ++col)
+            PPM_ASSIGN(pixels[leftR1][col], BG_RED, BG_GREEN, BG_BLUE);
+        for (col = 0; col < leftC2; ++col)
+            PPM_ASSIGN(pixels[leftR2][col], BG_RED, BG_GREEN, BG_BLUE);
 
-    procTop(c1, cm, r1, rm, var);
-    procTop(cm, c2, rm, r2, var);
+        procLeft(pixels, leftR1, leftR2, leftC1, leftC2, var, randStP);
+    }
+}
+
+
+
+static void
+makeRaggedRightBorder(pixel **           const pixels,
+                      unsigned int       const rows,
+                      unsigned int       const cols,
+                      unsigned int       const right,
+                      unsigned int       const width,
+                      unsigned int       const var,
+                      struct pm_randSt * const randStP) {
+
+    if (right >= 0) {
+        int const rightC1 = cols - right - 1;
+        int const rightC2 = cols - right - 1;
+        int const rightR1 = 0;
+        int const rightR2 = rows - 1;
+
+        unsigned int col;
+
+        for (col = rightC1; col < cols; ++col)
+            PPM_ASSIGN(pixels[rightR1][col], BG_RED, BG_GREEN, BG_BLUE);
+        for (col = rightC2; col < cols; ++col)
+            PPM_ASSIGN(pixels[rightR2][col], BG_RED, BG_GREEN, BG_BLUE);
+
+        procRight(pixels, rightR1, rightR2, rightC1, rightC2, width, var,
+                  randStP);
+    }
 }
 
 
 
 static void
-procBottom(int          const c1,
-           int          const c2,
-           int          const r1,
-           int          const r2,
-           unsigned int const height,
-           unsigned int const var) {
+makeRaggedTopBorder(pixel **           const pixels,
+                    unsigned int       const rows,
+                    unsigned int       const cols,
+                    unsigned int       const top,
+                    unsigned int       const var,
+                    struct pm_randSt * const randStP) {
 
-    int rm, cm, r;
+    if (top >= 0) {
+        unsigned int const topR1 = top;
+        unsigned int const topR2 = top;
+        unsigned int const topC1 = 0;
+        unsigned int const topC2 = cols - 1;
 
-    if (c1 + 1 == c2)  return;
-    cm = (c1 + c2) >> 1;
-    rm = (r1 + r2) >> 1;
-    rm += (int)floor(((float)rand() / RAND_MAX - 0.5) * var + 0.5);
+        unsigned int row;
 
-    for (r = rm; r < height; r++)
-        PPM_ASSIGN(PIX[r][cm], BG_RED, BG_GREEN, BG_BLUE);
+        for (row = 0; row < topR1; ++row)
+            PPM_ASSIGN(pixels[row][topC1], BG_RED, BG_GREEN, BG_BLUE);
+        for (row = 0; row < topR2; ++row)
+            PPM_ASSIGN(pixels[row][topC2], BG_RED, BG_GREEN, BG_BLUE);
 
-    procBottom(c1, cm, r1, rm, height, var);
-    procBottom(cm, c2, rm, r2, height, var);
+        procTop(pixels, topC1, topC2, topR1, topR2, var, randStP);
+    }
+}
+
+
+
+static void
+makeRaggedBottomBorder(pixel **            const pixels,
+                       unsigned int        const rows,
+                       unsigned int        const cols,
+                       unsigned int        const bottom,
+                       unsigned int        const height,
+                       unsigned int        const var,
+                       struct pm_randSt *  const randStP) {
+
+    if (bottom >= 0) {
+        unsigned int const bottomR1 = rows - bottom - 1;
+        unsigned int const bottomR2 = rows - bottom - 1;
+        unsigned int const bottomC1 = 0;
+        unsigned int const bottomC2 = cols - 1;
+
+        unsigned int row;
+
+        for (row = bottomR1; row < rows; ++row)
+            PPM_ASSIGN(pixels[row][bottomC1], BG_RED, BG_GREEN, BG_BLUE);
+        for (row = bottomR2; row < rows; ++row)
+            PPM_ASSIGN(pixels[row][bottomC2], BG_RED, BG_GREEN, BG_BLUE);
+
+        procBottom(pixels, bottomC1, bottomC2, bottomR1, bottomR2,
+                   height, var, randStP);
+    }
 }
 
 
@@ -188,14 +332,18 @@ main(int argc, const char * argv[]) {
     struct CmdlineInfo cmdline;
     pixel bgcolor, fgcolor;
     pixval fg_red, fg_green, fg_blue;
-    int rows, cols, row;
+    int rows, cols;
     int left, right, top, bottom;
+    struct pm_randSt randSt;
+    static pixel** pixels;
 
     pm_proginit(&argc, argv);
 
     parseCommandLine(argc, argv, &cmdline);
 
-    srand(cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
+    pm_randinit(&randSt);
+    pm_srand(&randSt,
+             cmdline.randomseedSpec ? cmdline.randomseed : pm_randseed());
 
     cols = cmdline.width;
     rows = cmdline.height;
@@ -236,92 +384,30 @@ main(int argc, const char * argv[]) {
         pm_message("foreground is %s",
                    ppm_colorname(&fgcolor, PPM_MAXMAXVAL, 1));
         if (cmdline.randomseedSpec)
-            pm_message("srand() initialized with seed %u", cmdline.randomseed);
-    }
-
-    /* Allocate memory for the whole pixmap */
-    PIX = ppm_allocarray(cols, rows);
-
-    /* First, set all pixel to foreground color */
-    for (row = 0; row < rows; row++) {
-        unsigned int col;
-        for (col = 0; col < cols; ++col)
-            PPM_ASSIGN(PIX[row][col], fg_red, fg_green, fg_blue);
-    }
-    /* Make a ragged left border */
-    if (left >= 0) {
-        int const left_c1 = left;
-        int const left_c2 = left;
-        int const left_r1 = 0;
-        int const left_r2 = rows - 1;
-
-        unsigned int col;
-
-        for (col = 0; col < left_c1; ++col)
-            PPM_ASSIGN(PIX[left_r1][col], BG_RED, BG_GREEN, BG_BLUE);
-        for (col = 0; col < left_c2; ++col)
-            PPM_ASSIGN(PIX[left_r2][col], BG_RED, BG_GREEN, BG_BLUE);
-
-        procLeft(left_r1, left_r2, left_c1, left_c2, cmdline.var);
-    }
-
-    /* Make a ragged right border */
-    if (right >= 0) {
-        int const right_c1 = cols - right - 1;
-        int const right_c2 = cols - right - 1;
-        int const right_r1 = 0;
-        int const right_r2 = rows - 1;
-
-        unsigned int col;
-
-        for (col = right_c1; col < cols; col++)
-            PPM_ASSIGN(PIX[right_r1][col], BG_RED, BG_GREEN, BG_BLUE);
-        for (col = right_c2; col < cols; col++)
-            PPM_ASSIGN(PIX[right_r2][col], BG_RED, BG_GREEN, BG_BLUE);
-
-        procRight(right_r1, right_r2, right_c1, right_c2,
-                   cmdline.width, cmdline.var);
+            pm_message("pm_rand() initialized with seed %u",
+                       cmdline.randomseed);
     }
 
-    /* Make a ragged top border */
-    if (top >= 0) {
-        int const top_r1 = top;
-        int const top_r2 = top;
-        int const top_c1 = 0;
-        int const top_c2 = cols - 1;
+    pixels = ppm_allocarray(cols, rows);
 
-        unsigned int row;
+    makeAllForegroundColor(pixels, rows, cols, fg_red, fg_green, fg_blue);
 
-        for (row = 0; row < top_r1; ++row)
-            PPM_ASSIGN(PIX[row][top_c1], BG_RED, BG_GREEN, BG_BLUE);
-        for (row = 0; row < top_r2; ++row)
-            PPM_ASSIGN(PIX[row][top_c2], BG_RED, BG_GREEN, BG_BLUE);
+    makeRaggedLeftBorder(pixels, rows, cols, left, cmdline.var, &randSt);
 
-        procTop(top_c1, top_c2, top_r1, top_r2, cmdline.var);
-    }
+    makeRaggedRightBorder(pixels, rows, cols, right,
+                          cmdline.width, cmdline.var, &randSt);
 
-    /* Make a ragged bottom border */
-    if (bottom >= 0) {
-        int const bottom_r1 = rows - bottom - 1;
-        int const bottom_r2 = rows - bottom - 1;
-        int const bottom_c1 = 0;
-        int const bottom_c2 = cols - 1;
+    makeRaggedTopBorder(pixels, rows, cols, top, cmdline.var, &randSt);
 
-        unsigned int row;
-
-        for (row = bottom_r1; row < rows; ++row)
-            PPM_ASSIGN(PIX[row][bottom_c1], BG_RED, BG_GREEN, BG_BLUE);
-        for (row = bottom_r2; row < rows; ++row)
-            PPM_ASSIGN(PIX[row][bottom_c2], BG_RED, BG_GREEN, BG_BLUE);
+    makeRaggedBottomBorder(pixels, rows, cols, bottom,
+                           cmdline.height, cmdline.var, &randSt);
 
-        procBottom(bottom_c1, bottom_c2, bottom_r1, bottom_r2,
-                   cmdline.height, cmdline.var);
-    }
+    pm_randterm(&randSt);
 
     /* Write pixmap */
-    ppm_writeppm(stdout, PIX, cols, rows, PPM_MAXMAXVAL, 0);
+    ppm_writeppm(stdout, pixels, cols, rows, PPM_MAXMAXVAL, 0);
 
-    ppm_freearray(PIX, rows);
+    ppm_freearray(pixels, rows);
 
     pm_close(stdout);
 
@@ -329,4 +415,3 @@ main(int argc, const char * argv[]) {
 }
 
 
-
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..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
+};
+
+