Mps br/timing rnd():

Use a = 48271.  Use fast implementation.  Much improved verification by rnd_verify().

Copied from Perforce
 Change: 167182
 ServerID: perforce.ravenbrook.com

2 files changed, 155 insertions, 54 deletions

diff --git a/mps/code/testlib.c b/mps/code/testlib.c
index 5c92969d2e0..947ead5fb82 100644
--- a/mps/code/testlib.c
+++ b/mps/code/testlib.c
@@ -34,72 +34,167 @@ struct itimerspec; /* stop complaints from time.h */
 
 /* rnd -- a random number generator
  *
- * I nabbed it from "ML for the Working Programmer", originally from:
- * Stephen K Park & Keith W Miller (1988). Random number generators:
- * good ones are hard to find.  Communications of the ACM, 31:1192-1201.
+ * We use the (Multiplicative) Linear Congruential Generator
+ *   Xn = a * Xn-1 mod m
+ * with: m = 2147483647 (2^31 - 1, a Mersenne prime), and a = 48271.  
+ * This is a 'full-period' generator: all values from [1..(mod-1)], 
+ * or 0x00000001 to 0x7ffffffe inclusive, are returned once, and then 
+ * the cycle begins again.  (The value 0 is not part of the cycle and 
+ * is never returned).  So the period = mod-1, ie. 2147483646.
  *
- * This is a 'full-period' generator: all values from [1..(mod-1)] 
- * are returned once, and then the generator cycles round again.
- * (The value 0 is not part of the cycle and is never returned). That 
- * means the period = mod-1, ie. 2147483646.
+ * This generator is extremely simple and has been very well studied.  
+ * It is free of major vices we might care about for this application.
+ * In particular, as m is prime, low order bits are random.  Therefore 
+ * to roll an N-sided die (N << m), "rnd() % N" is acceptable, giving 
+ * a value in [0..N-1].
  *
- * This "x 16807 mod 2^31-1" generator has been very well studied.  
- * It is not the 'best' (most random) simple generator, but it is 
- * free of all vices we might care about for this application, so 
- * we'll keep it simple and stick with this, thanks.
+ * It was popularised by the much-cited Park & Miller paper:
+ *   Stephen K Park & Keith W Miller (1988). Random number generators:
+ *   good ones are hard to find.  Communications of the ACM, 
+ *   31:1192-1201.
+ * The recommended multiplier a was later updated from 16807 to 48271:
+ *   Stephen K Park, Keith W Miller, Paul K. Stockmeyer (1993). 
+ *   Technical Correspondence.  Communications of the ACM, 36:105-110.
  *
- * There are faster implementations of this generator, summarised 
- * briefly here:
- *  - L. Schrage (1979 & 1983) showed how to do all the calculations 
- *    within 32-bit integers.  See also Numerical recipes in C.
- *  - David Carta (1990) showed how to also eliminate division.  
+ * (Many more elaborate generators have been invented.  The next simple 
+ * step would be to combine with the MLCG m = 2147483399 a = 40692, to 
+ * make the period "about 74 quadrillion".  See the summary of chapter 
+ * 3 in Knuth's "The Art of Computer Programming".)
+ *
+ * This (fast) implementation uses the identity:
+ *   0x80000000 == 0x7FFFFFFF + 0x00000001
+ * noted by David Carta (1990), where 0x7FFFFFFF == 2^31-1 == m, which 
+ * means that bits above the first 31 can simply be shifted >> 31 and 
+ * added, preserving Xn mod m.  To remain within 32-bit unsigned 
+ * arithmetic when multiplying the previous seed (31 bits) by a (16 
+ * bits), the seed is split into bottom and top halves; bits above 
+ * the first 31 are simply "top >> 16".  (Code by RHSK, inspired by 
+ * Robin Whittle's article at "http://www.firstpr.com.au/dsp/rand31/").
+ *
+ * Slower implementations, used for verification:
+ * rnd_verify_schrage uses the method of L. Schrage (1979 & 1983), 
+ *   namely splitting the seed by q, where q = m div a.
+ * rnd_verify_float simply uses floating point arithmetic.
  */
 
 static unsigned long seed = 1;
-/* 2^31 - 1 == (1UL<<31) - 1 == 2147483647 */
-#define RND_MOD_2_31_1       2147483647
-#define RND_MOD_2_31_1_FLOAT 2147483647.0
+#define R_m 2147483647UL
+#define R_a 48271UL
 unsigned long rnd(void)
 {
-  double s;
-
-  s = seed;
-  s *= 16807.0;
-  s = fmod(s, RND_MOD_2_31_1_FLOAT);  /* 2^31 - 1 */
-  seed = (unsigned long)s;
+  /* requires m == 2^31-1, a < 2^16 */
+  unsigned long bot = R_a * (seed & 0x7FFF);
+  unsigned long top = R_a * (seed >> 15);
+  seed = bot + ((top & 0xFFFF) << 15) + (top >> 16);
+  if(seed > R_m)
+    seed -= R_m;
   return seed;
-  /* Have you modified this code?  Run rnd_verify() please!  RHSK */
+  /* Have you modified this code?  Run rnd_verify(3) please!  RHSK */
 }
 
-/* rnd_verify -- verify that rnd() returns the correct results */
-static void rnd_verify(void)
+static unsigned long seed_verify_schrage = 1;
+#define R_q (R_m / R_a)
+#define R_r (R_m % R_a)
+static unsigned long rnd_verify_schrage(void)
+{
+  /* requires m < 2^31, q > r; see Park & Miller (1988) */
+  unsigned long alpha = R_a * (seed_verify_schrage % R_q);  /* < m */
+  unsigned long beta = R_r * (seed_verify_schrage / R_q);  /* < m */
+  seed_verify_schrage = alpha - beta;
+  if(alpha < beta)
+    seed_verify_schrage += R_m;
+  return seed_verify_schrage;
+}
+
+static unsigned long seed_verify_float = 1;
+#define R_m_float 2147483647.0
+#define R_a_float 48271.0
+static unsigned long rnd_verify_float(void)
+{
+  double s;
+  s = seed_verify_float;
+  s *= R_a_float;
+  s = fmod(s, R_m_float);
+  seed_verify_float = (unsigned long)s;
+  return seed_verify_float;
+}
+
+/* rnd_verify -- verify that rnd() returns the correct results
+ *
+ * depth = how much time to spend verifying
+ * 0: very quick -- just verify the next rnd() value
+ * 1: quick -- verify the first 10000 calls from seed = 1
+ * 2: slow (< 5 minutes) -- run the fast generator for a full cycle
+ * 3: very slow (several hours) -- verify a full cycle
+ */
+void rnd_verify(int depth)
 {
   unsigned long orig_seed = seed;
   unsigned long i;
+  unsigned long r;
+  
+  /* 0: the next value from rnd() matches rnd_verify_*() */
+  if(depth >= 0) {
+    seed_verify_schrage = seed;
+    seed_verify_float = seed;
+    r = rnd();
+    Insist(rnd_verify_schrage() == r);
+    Insist(rnd_verify_float() == r);
+  }
 
-  /* This test is from:
-   * Stephen K Park & Keith W Miller (1988). Random number generators:
-   * good ones are hard to find.  Communications of the ACM, 31:1192-1201.
-   * Note: The Park-Miller paper uses 1-based indices.
-   */
-  i = 1;
-  seed = 1;
-  for(i = 2; i < 10001; i += 1) {
-    (void)rnd();
+  /* 1: first 10000 (from Park & Miller, note: 1-based indexing!) */
+  if(depth >= 1) {
+    i = 1;
+    seed = 1;
+    seed_verify_schrage = 1;
+    seed_verify_float = 1;
+    for(i = 2; i <= 10001; i += 1) {
+      r = rnd();
+      Insist(rnd_verify_schrage() == r);
+      Insist(rnd_verify_float() == r);
+    }
+    printf("%lu\n", r);
+    /* Insist(r == 1043618065UL); -- correct value for a = 16807 */
+    Insist(r == 399268537UL);  /* correct for a = 48271 */
   }
-  Insist(i == 10001);
-  Insist(rnd() == 1043618065);
 
-  /* from Robin Whittle's compendium at 
-   * http://www.firstpr.com.au/dsp/rand31/
-   * Observe wrap-around after 'full-period' (which excludes 0).
-   * Note: uses 0-based indices, ie. rnd_0 = 1, rnd_1 = 16807, ...
-   */
-  i = 2147483645;
-  seed = 1407677000;
-  i += 1;
-  Insist(i == (1UL<<31) - 2 ); /* 'full-period' excludes 0 */
-  Insist(rnd() == 1);  /* wrap-around */
+  /* 1: observe wrap-around (note: 0-based indexing) */
+  if(depth >= 1) {
+    i = 2147483645UL;
+    /* seed = 1407677000UL; -- correct value for a = 16807 */
+    seed = 1899818559UL;  /* correct for a = 48271 */
+    seed_verify_schrage = 1899818559UL;
+    seed_verify_float = 1899818559UL;
+    i += 1;
+    Insist(i == (1UL<<31) - 2 ); /* 'full-period' excludes 0 */
+    Insist(rnd() == 1);  /* wrap-around */
+    Insist(rnd_verify_schrage() == 1);
+    Insist(rnd_verify_float() == 1);
+  }
+
+  /* 2 & 3: Full cycle (3 => verifying each value) */
+  if(depth >= 2) {
+    unsigned long r1;
+    i = 0;
+    seed = 1;
+    seed_verify_schrage = 1;
+    seed_verify_float = 1;
+    while(1) {
+      i += 1;
+      r = rnd();
+      if(depth >= 3) {
+        Insist(rnd_verify_schrage() == r);
+        Insist(rnd_verify_float() == r);
+      }
+      if(r == 1) {
+        printf("Wrap at i=%lu, r=%lu, r(i-1)=%lu.\n",
+               i, r, r1);
+        break;
+      } else {
+        r1 = r;
+      }
+    }
+  }
 
   seed = orig_seed;
 }
@@ -135,18 +230,20 @@ void randomize(int argc, char **argv)
     n = sscanf(argv[1], "%lu", &seed0);
     Insist(n == 1);
     printf("randomize(): resetting initial seed to: %lu.\n", seed0);
-    rnd_verify();  /* good to verify occasionally: slip it in here */
   } else {
     /* time_t uses an arbitrary encoding, but hopefully the low order */
     /* 31 bits will have at least one bit changed from run to run. */
-    seed0 = 1 + time(NULL) % (RND_MOD_2_31_1 - 1);
+    seed0 = 1 + time(NULL) % (R_m - 1);
     printf("randomize(): choosing initial seed: %lu.\n", seed0);
   }
 
-  Insist(seed0 < RND_MOD_2_31_1);  /* 2^31 - 1 */
+  Insist(seed0 < R_m);
   Insist(seed0 != 0);
   seed = seed0;
 
+  rnd_verify(0);
+  Insist(seed == seed0);
+
   /* The 'random' seed is taken from time(), which may simply be a 
    * count of seconds: therefore successive runs may start with 
    * nearby seeds, possibly differing only by 1.  So the first value
diff --git a/mps/code/testlib.h b/mps/code/testlib.h
index 94161adb58d..cf40053c2d3 100644
--- a/mps/code/testlib.h
+++ b/mps/code/testlib.h
@@ -138,12 +138,16 @@ extern void verror(const char *format, va_list args);
 
 /* rnd -- random number generator
  *
- * rnd() generates a sequence of integers in the range [0, 2^31-2].
+ * rnd() generates a sequence of integers in the range [1, 2^31-1].
  */
 
 extern unsigned long rnd(void);
 
 
+/* rnd_verify() -- checks behaviour of rnd() */
+extern void rnd_verify(int depth);
+
+
 /* rnd_addr -- random number generator
  *
  * rnd_addr() generates a sequence of addresses all over the address space.
@@ -167,7 +171,7 @@ extern void randomize(int argc, char **argv);
 
 /* C. COPYRIGHT AND LICENSE
  *
- * Copyright (C) 2001-2002 Ravenbrook Limited <http://www.ravenbrook.com/>.
+ * Copyright (C) 2001-2002, 2008 Ravenbrook Limited <http://www.ravenbrook.com/>.
  * All rights reserved.  This is an open source license.  Contact
  * Ravenbrook for commercial licensing options.
  *