From e486078a79edd55a290e4087ee92c114fa301e0d Mon Sep 17 00:00:00 2001 From: Richard Kistruck Date: Tue, 20 Jan 2009 15:59:32 +0000 Subject: [PATCH] Mps br.timing: testlib randomize() was wrong (simply discarding the first 0..32000 values is a throughly defective way to seed the generator). correct it to simply use timestamp (modulo m) as the initial seed, and then iterate 10 times. rnd(): add notes, and create rnd_verify() to confirm implementation is correct. Copied from Perforce Change: 167181 ServerID: perforce.ravenbrook.com --- mps/code/testlib.c | 99 ++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 87 insertions(+), 12 deletions(-) diff --git a/mps/code/testlib.c b/mps/code/testlib.c index b22ed52bd2c..5c92969d2e0 100644 --- a/mps/code/testlib.c +++ b/mps/code/testlib.c @@ -36,21 +36,73 @@ struct itimerspec; /* stop complaints from time.h */ * * I nabbed it from "ML for the Working Programmer", originally from: * Stephen K Park & Keith W Miller (1988). Random number generators: - * good one are to find. Communications of the ACM, 31:1192-1201. + * good ones are hard to find. Communications of the ACM, 31:1192-1201. + * + * 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 "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. + * + * 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. */ +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 unsigned long rnd(void) { - static unsigned long seed = 1; double s; s = seed; s *= 16807.0; - s = fmod(s, 2147483647.0); /* 2^31 - 1 */ + s = fmod(s, RND_MOD_2_31_1_FLOAT); /* 2^31 - 1 */ seed = (unsigned long)s; return seed; + /* Have you modified this code? Run rnd_verify() please! RHSK */ } +/* rnd_verify -- verify that rnd() returns the correct results */ +static void rnd_verify(void) +{ + unsigned long orig_seed = seed; + unsigned long i; + + /* 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(); + } + 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 */ + + seed = orig_seed; +} /* rnd_addr -- a random address generator * @@ -75,19 +127,42 @@ mps_addr_t rnd_addr(void) void randomize(int argc, char **argv) { - int i, k, n; + int n; + unsigned long seed0; + int i; if (argc > 1) { - n = sscanf(argv[1], "%d", &k); - die((n == 1) ? MPS_RES_OK : MPS_RES_FAIL, "randomize"); + 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 { - k = time(NULL) % 32000; - printf("Randomizing %d times.\n", k); + /* 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); + printf("randomize(): choosing initial seed: %lu.\n", seed0); } - /* Randomize the random number generator a bit. */ - for (i = k; i > 0; --i) - rnd(); + Insist(seed0 < RND_MOD_2_31_1); /* 2^31 - 1 */ + Insist(seed0 != 0); + 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 + * returned by rnd() may differ by only 16807. It is conceivable + * that some tests might be able to 'spot' this pattern (for + * example: by using the first rnd() value, mod 100M and rounded + * to multiple of 128K, as arena size). + * + * So to mix it up a bit, we do a few iterations now. How many? + * Very roughly, 16807^2 is of the same order as 2^31, so two + * iterations would make the characteristic difference similar to + * the period. Hey, let's go wild and do 10. + */ + for(i = 0; i < 10; i += 1) { + (void)rnd(); + } } @@ -146,7 +221,7 @@ void cdie(int res, const char *s) /* C. COPYRIGHT AND LICENSE * - * Copyright (C) 2001-2002 Ravenbrook Limited . + * Copyright (C) 2001-2002, 2008 Ravenbrook Limited . * All rights reserved. This is an open source license. Contact * Ravenbrook for commercial licensing options. *