Use Mersenee Twister instead of libc rand for random numbers

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

6 files changed, 676 insertions, 452 deletions

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[]) {
 }
 
 
-