/* newrand32.c 0.03 dh:2012-12-02 * * DEMONSTRATION OF 32-BIT WICHMANN-HILL PSEUDO-RANDOM GENERATOR * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ #include #include #include const long modx = 2147483579L; const long mody = 2147483543L; const long modz = 2147483423L; const long modt = 2147483123L; static long ix = 2080630037L, iy = 2033578565L, iz = 2101153635L, it = 2129003085L; /* initializing state in case there's no call on snewrandom32() first. */ double newrand32( ) { /* This is the updated Wichman-Hill PRG with near-31-bit precision. The implementation is based on the variant described in [Wichmann, B.A., Hill, I.D. Generating good pseudo-random numbers. Computational Statistics & Data Analysis 51, 3 (December 2006), pp. 1614-1622] Available at ]. This code is based on that in the December 5, 2005, long- version draft available from eurometros.org at . The technique for computing (a*x) mod n without requiring 64-bit arithmetic is also presented in Problems 3.2.1.1(9)-(11) and their solutions in [Knuth, Donald E. The Art of Computer Programming, vol.2: Seminumerical Algorithms, ed.3. Addison-Wesley (Reading MA: 1998)]. */ long double W; ix = 11600L*(ix%185127L) - 10379L*(ix/185127L); if (ix < 0) ix += modx; iy = 47003L*(iy%45688L) - 10479L*(iy/45688L); if (iy < 0) iy += mody; iz = 23000L*(iz%93368L) - 19423L*(iz/93368L); if (iz < 0) iz += modz; it = 33000L*(it%65075L) - 8123L*(it/65075L); if (it < 0) it += modt; W = (long double)ix / modx + (long double)iy / mody + (long double)iz / modz + (long double)it / modt; return (double)(W - (long)W); /* XXX: Assuming that casting down from long double will never cause 1.0 to be produced. */ } /* newrand32 */ static long simpleR(unsigned long *seeder, long whmod) { /* Simple PRG used for expansion of a single seed value to four initial seeds for the Whitmann-Hill 32-bit PRG. Generate pseudo-randoms, R, in the required range 0 < R < whmod < 2^31 and whmod > 2^30. This uses the Marsaglia "best possible multiplier candidate" simple modular generator (69069 * x) mod 2^32. This is the default seed generator for Mersenne Twister adapted here to limit the range of R to less than whmod. See also Table 1 in section 3.3.4 of [Knuth, Donald E. The Art of Computer Programming, vol.2: Seminumerical Algorithms, ed.3, Addison-Wesley (Reading, MA: 1998)]. */ unsigned long r; r = *seeder; while(1) { r = (unsigned long)((r * 69069uL) & 0xFFFFFFFFuL); if (r < (unsigned long)whmod && r) return (*seeder = r); } } /* simpleR */ void snewrand32(unsigned long seed) {/* Given a sufficiently unpredictable seed, generate the initial seeding of the Whitmann-Hill 32-bit pseudo-random generator. */ int scrambler = 1023; /* Make large enough that it is costly to find a seed that provokes degenerate results from the first calls on newrand32() */ unsigned long R; do { R = ((seed ? seed : 0xA5A55A5AuL) * 0x5A5AA5A5uL) & 0xFFFFFFFFuL;} while (!R); /* Scramble the seed to spread apart similar values. Prevent pathological 0 value for R. Already-random seeds will work just fine. See the discussion around condition 6.4(6) in [Knuth, Donald E. The Art of Computer Programming, vol.3: Sorting and Searching, ed.2, Addison-Wesley (Reading, MA: 1998)]. */ iy = simpleR(&R,mody); ix = simpleR(&R,modx); it = simpleR(&R,modt); iz = simpleR(&R,modz); /* Re-ordering just to be different from other literal implementations. */ while (scrambler--) (void) newrand32(); /* Counteracting degenerate values of R that degrade initial results. */ } /* snewrand32 */ int main(int argc, char **argv) { long i = 0; long run = 5; double r; if (argc > 1) run = atol(argv[1]); snewrand32( (unsigned long)(time(NULL)) ); /* XXX: This is not a great seed. It is used here to show that the scrambler in snewrand32 is effective regardless. */ for (i = 0; i < run; i++) { r = newrand32(); printf("%23.20f%12ld:x%12ld:y%12ld:z%12ld:t\n", r, ix, iy, iz, it); } return 0; } /* main */ /* 0.03 2012-12-03T02:10Z Update snewrandom32() to block pathological cases of seeds that scramble to 0 or to trivial values that produce unsatisfactory initial newrandom32() results. It's all heuristic. 0.02 2012-12-02T20:49Z Replaced the clumsy scrambling of time(NULL) in the call to snewrand32() with a good scrambler in snewrand32() that should work well for all qualities of user-provided seeds. 0.01 2012-12-02T04:36Z Uploaded version with snewrandom32( ) and a simpleR technique for bootstrapping initial seeds via a simple modular 32-bit generator. 0.00 2012-12-01T18:36Z Created a 31-bit version based on the 2005 improvement of the Wichmann-Hill method using 4 cycles that is far superior to the original 15-bit version. This was inspired by back-and-forth discussions and trial implementations with Pedro Guffini while he was looking to provide a better RAND() function for Apache OpenOffice.calc. */ /* end of newrand32.c */