Abstract: Pseudo-random number generator by Warren D. Smith.
Goal is to combine several random number generators to get at least the randomness of the strongest of them, but to do it in such a way as to get nearly the speed of the fastest of them.
NOTE: I use a lot of global variables here to save state,
necessary in modula3 since C's
static feature does not exist.
Of course these variables are local to the current module.
For that reason, please keep this a separate module from others,
e.g. the derived routines for NormalDeviate, etc. should not be in
same module as these underlying routines for uniform deviates.
Some component generators: Following G.Marsaglia: A current view of RNGs, pages 3-10 in Computer Science and Statistics, the interface, Elsevier 1985, a fair number of my generators are of the FIBO(a,b,%) type, where a>b>0 are integers such that x^a+x^b+1 is a primitive trinomial mod 2 and % is a binary operation.
Suitable values of (b,a) are tabulated in D.Knuth: Seminumerical methods,
Addison-Wesley 1981, page 28, which is extracted from larger tables in
N.Zierler & J.Brillhart: Information and Control 13 (1968) 541-554;
14 (1969) 566-569; 15 (1969) 67-69. Some pairs with
a prime which
go beyond the Knuth table, which you can use if you *really*
want randomness, include:
(7,127), (15,127), (30,127), (63,127), (32,521), (158,521),
(105,607), (147,607), (273,607), (216,1279), (418,1279)
The binary operation % could be: [- mod 1] for floating point values, [- mod M] or [* mod M] for some integer M, or [this NOT recommended; it is same as - only mod 2] XOR.
Then the meta-procedure Fibo(a,b,%) is as follows (it requires an auxiliary arr : ARRAY [0..a-1], initially filled with random seed values, and initially ia=a and ja=b:
MetaRandGen() : Type = BEGIN DEC(i); DEC(j); IF i<0 THEN i := a-1; (* wraparound
ELSIF j<0 THEN j := a-1; (* wraparound *) END; arr[i] := arr[i] % arr[j]; RETURN arr[i]; END MetaRandGen; Period [Marsaglia & Tsay: Lin. Alg. and its Appl. 67, 147-155, 1985]: If % is (- mod 2^wordsize), then if at least one of the arr is odd the period will be maximal = (2^a - 1) * 2^(wordsize-1) if the characteristic axa matrix T of the recurrence has full order J=(2^a - 1) in mod 2 arithmetic (this assured by the trinomial condition) and 2*J in mod 4 and 4*J in mod 8 arithmetic. That is, if you square T a times, you get back T in mod 2, but do NOT get back T in mod 4, arithmetic, and if you square T a+1 times, you do not get back T in mod 8 arithmetic. In fact, this criterion will work for any matrix T, not necessarily the fibo(-) one... If % is [* mod 2^wordsize], which may be implemented in m3 by Word.Times(), then arr must be all odd integers, and by considering discrete logs we see that the period is maximal = (2^a - 1) * 2^(wordsize-3) if the F(a,b,+ mod 2^wordsize) has maximal period (2^a - 1) * 2^(wordsize-1). Marsaglia's test results for Fibo(a,b,%) generators with (a,b)=(17,5), (31,13) and (55,24): XOR: With %=XOR the fibo fail seven of the tests in Marsaglia's battery. As Marsaglia says in his conclusion "never use XOR." Despite this, these generators keep resurfacing, for example the "r250" Fibo(250,147,XOR) generator of S. Kirkpatrick and E. Stoll, Journal of Computational Physics, 40, p. 517 (1981) and W. L. Maier, "A Fast Pseudo Random Number Generator", Dr. Dobb's Journal #176. Subtraction: With %=-, among the tests in Marsaglia's DIEHARD battery only the "birthday spacings test" (of frequencies of spacings in sorted sets of deviates) is failed, a symptom both of the linear structure and also of the specific subtractive structure of these generators. Note that the Fibo(55,24,-) generator (which fails this test) is to be found in Knuth's book and also is distributed with the DEC SRC implementation of modula-3. Also note that postprocessing Fibo(55,24,-) sequence with a Bays-Durham shuffler, recommended in Press et al. "Numerical Recipes" to fix suspicious generators, will NOT work to make it now pass this test, since the birthday spacings test does not depend upon the ordering of the deviates only their values. Multiplication: With %=* mod 2^32, these FIBO gens passed all the tests in Marsaglia's battery. Note that the -,+ and XOR gens are "linear" and hence theoretically bad and will always fail "empty slab tests" see below. Note that shift register, Fibonacci() generators with +,- or XOR, and linear congruential generators x <-- a*x+b mod M, AND linear combinations of such possibly to different moduli, are all "linear" and thus generate d-tuples of random numbers lying on AT MOST M^(1/d) hyperplanes in d-space. Consequently, the nonrandomness of a linear RNG with period M is in principle detectable by the "largest empty slab" test in .4*logM dimensions after only N, N > e^(5/2) * (.4*logM)^2.5 * loglogM, random numbers have been generated. Hence: Please do not rely on a linear generator. If you are going, foolishly, to use a linear congruential generator, though, you want a "good multiplier" a mod M. (Bad multipliers will result in even fewer hyperplanes, for example IBM's infamous RANDU generator generated points lying on only 15 planes in 3-space. Even with fairly good multipliers we still get test failures, e.g. Marsaglia's spectrally good multiplier 69069 mod 2^32 fails his OPSO test, as does the Berkeley PASCAL 62605 mod 2^29.) It will suffice, for comparatively good behavior in d-space, that (a^d mod M)/M have only small partial quotients in its continued fraction expansion. This should be tested for d=1,2,..,8 at least. Thus for example, consider the prime modulus M=2^35 - 849. (M-1)/2 is also prime. I found the generator g=145683 by computer search. The continued fraction expansions CF(g/p)=[235852,1,3,6,1,1,5,1,1,2,1,1,1,1,6,2] CF(g^2/p)=[1,1,1,1,1,1,1,1,21,2,7,1,8,1,2,1,1,2,5,1,1,1,2,1,2,3,12,2] are rather nice, and the CF( (g^x mod p)/p ) for x=3,4,5,6,7,8,9 are also not bad (no partial quotient larger than 23): [2,1,1,1,1,1,1,1,1,16,1,3,5,1,2,1,1,3,1,1,1,21,1,8,10,5] [8,4,3,1,1,3,4,2,7,1,22,22,9,1,12,1,2,2] [1,1,1,20,1,23,1,4,2,2,1,1,1,8,1,1,2,1,2,2,1,22,1,2,1,3] [1,1,15,1,1,10,2,2,2,1,4,5,1,2,1,1,5,7,7,2,2,4,2] [1,1,5,2,2,1,1,1,7,11,3,4,1,9,1,6,3,1,8,1,3,1,1,1,2] [1,13,1,7,8,1,1,2,16,2,2,1,1,6,4,1,3,1,1,3,15,2] [2,1,2,1,15,2,1,5,1,8,1,1,2,1,4,3,1,3,1,1,1,1,2,1,3,1,3,5] so we conclude that using g as a multiplier should exhibit comparatively good behavior in dimensions 2-9. This particular (g,M) pair has the virtue that IEEE doubles can represent integers up to and including 2^53 - 1 exactly, so that the modula-3 statement x := (g*x) MOD M; will evaluate it exactly. As another example, the Marsaglia multiplier 69069 mod 2^32, while spectrally good in dimensions 2-5, is bad in dimension 6, as is revealed by the spectral test directly but also simply by noticing that the CF expansion of 69069^6 / 2^32 is [1, 75, 1, 2, 2, 1, 20, 10, 3, 10, 2, 2, 1, 12, 9] which contains the large number 75 early on. You also probably want full period M, which happens if [thm page 16 Knuth] (1) GCD(b,M)=1; (2) a-1 is a multiple of p for every prime p dividing M; (3) a-1 is a multiple of 4 if 4 divides M. If M is prime and b=0, you get maximal period M-1 if a primitive mod M. If M>=16 is power of 2 and b=0, get maximal period M/4 if a=3 or 5 mod 8. Other Fibo generators: instead of just using one lag term and binary operation %, you could combine with TWO lag terms via some TERNARY operation, or THREE lag terms via a QUATERNARY operation, etc. I suggest the new generator QuaternaryFibo(a,b,c,d, x0 - (x2 XOR (x1 - x3) mod M) mod M ). If M is 2^wordsize and x^a+x^b+x^c+x^d+1, a>b>c>d>0, is a primitive polynomial mod 2, then this generator's period will be at least 2^a - 1 simply by considering the LS bit, which will follow a DeBruijn sequence and thus exhibits good randomness in <=a dimensions. My idea is that by using both XOR and +, we hope to avoid the weak behavior of either operation alone, e.g. with respect to the birthday spacings test. This may also make the generator "nonlinear" (i.e. not outputting a lattice), a point I am unsure about. Ternary generators of this type do not seem to exist since 4-term primitive polynomials x^a+x^b+x^c+1 mod 2 do not exist (you need an odd number of terms, as can easily be shown). For tables of examples of primitive polynomials mod 2, see E.J.Watson: Math. of Comput. 16 (1962) 368-9 but the ones given are the lexically first examples, bad for our purposes. Hence I constructed my own small table below of examples of primitive polynomials P=x^a+x^b+x^c+x^d+1 mod 2, where a is prime and N=2^a - 1 is a Mersenne prime. For such a, P is primitive if P^N divides x^N-1 mod 2, which you can test by the algorithm Q := x; FOR i:=1 TO a-1 do Q := Q*Q*x mod P; END; IF Q=1 THEN "P is primitive." Table: x^19 + x^9 + x^4 + x^3 + 1 x^19 + x^6 + x^2 + x^1 + 1 x^31 + x^22 + x^20 + x^8 + 1 x^31 + x^16 + x^11 + x^6 + 1 x^61 + x^5 + x^2 + x^1 + 1 x^61 + x^57 + x^2 + x^1 + 1 x^61 + x^19 + x^16 + x^9 + 1 x^89 + x^36 + x^2 + x^1 + 1 x^89 + x^57 + x^14 + x^5 + 1 x^89 + x^28 + x^26 + x^9 + 1 x^89 + x^37 + x^5 + x^1 + 1 x^107 + x^53 + x^45 + x^12 + 1 x^107 + x^93 + x^68 + x^2 + 1 x^127 + x^29 + x^17 + x^7 + x^127 + x^67 + x^65 + x^17 + 1 x^127 + x^60 + x^22 + x^9 + 1 x^521 + x^46 + x^38 + x^37 + 1 x^521 + x^249 + x^92 + x^87 + 1 x^521 + x^353 + x^258 + x^6 + 1 x^607 + x^66 + x^60 + x^30 + 1 x^607 + x^247 + x^83 + x^71 + 1 x^607 + x^483 + x^449 + x^298 + 1 x^607 + x^442 + x^113 + x^23 + 1 It is a disgusting fact that in modula3 as well as most high level languages, the machine language multiply instruction that computes a*b to double precision is NOT ACCESSIBLE. Ditto "add with carry". This makes an efficient implementation of high precision multiplication (& modular version) obnoxiously difficult and inefficient. And you really do need multiple precision to get a decent period with a lincong or iterated squaring generator. [Note the SQUARE ROOT of the period must be unreachably large for good randomness, NOT the period, a common error.] Here is a MetaAlgorithm which works if numbers 0..modulus*2-2 are representable...: MetaModularMultiply(x,y, modulus: NumType) : NumType = VAR q : NumType := x; BEGIN FOR b = bits of y in MS-->LS order starting at the bit AFTER the first 1 bit DO q := q+q; IF q>=modulus THEN q := q-modulus; END; IF b=1 THEN q := q+x; IF q>=modulus THEN q := q-modulus; END; END; END; RETURN q; END MetaModularMultiply; One way to improve the randomness of a generator is C.Bays and S.D. Durham's shuffling algorithm, ACMTOMS 2 (1976) 59-64, also described in D.Knuth: Seminumerical algorithms. Knuth says this will output a "considerably more random" sequence. However, the set of values will still be the same (although ordered differently) hence, e.g., the subtractive and XOR fibo gens will STILL fail Marsaglia's birthday spacings test even after improvement by shuffling. However, linear congruential generators are probably substantially improved by shuffling since the lattice structure is destroyed. L.Blum, M.Blum, M. Shub: A simple unpredictable psuedo-random number generator, SIAM J. Comput 15,2 (1986) 364- following A.Yao: Theory and applications of trapdoor fns, 23rd SFOCS (1982) 80-91 and W.Alexi, B.Chor, O.Goldreich, C.P.Schnorr: 25th SFOCS (1984) 449-457. point out that iterated squaring generator x <--- x*x mod M where M=p*q, p, q primes that are 3 mod 4, will have these properties (assuming it is computationally infeasible to compute the factorization of M) (1) infeasible to predict the PREVIOUS x. (2) infeasible to predict the least signif bit of the previous x with correctness prob > 1/2 + 1/polynomial. (3) infeasible to predict the boolean "x<(M-1)/2" for the previous x with correctness prob > 1/2 + 1/polynomial. consequence of (2) is: no polynomial time statistical test can invalidate the randomness of the sequence of LS bits of the iterated squaring generator. The full iterated squaring gen, not just its LS bits, is a lot faster and (conjecturally!) just as random. Its period is (p-1)*(q-1)/4. Example: Here is a product of two primes both 3 mod 4: 94906247 * 94906219 = 9007193062250093 = 2^53 - 6192490899. Here is a smaller example: M = 9739 * 9719 = 94653341. This modulus has the virtue that iterated squaring mod M may be accomplished exactly in IEEE double arithmetic using x := x*x MOD M, since M^4 < 2^53 - 1. An even smaller modulus suitable for 32-bit unsigned arithmetic is M = 239*251; (* = 59989; 239 and 251 are each primes and 3 mod 4 *). More generally you could use iterated cubing (subject to the conjecture it is infeasible to break any bit of the RSA cryptosystem) or in fact any fixed exponent you want (under same conjecture) and in fact almost any fixed integer linear operation on the discrete logarithms of a vector will by the same reasoning be immune to any polytime statistical test, subject to RSA-type conjectures. In particular, we conjecture that the Fibo[a,b,*] generator ought to be immune to polytime statistical tests if the modulus is the product of two suitable large primes... which explains the results of the Marsaglia test battery above. 3/17/96 Warren Smith Initial version 3/23/96 Harry George Added object wrappers, which also required "InitDone" flag in Init proc. *) IMPORT LongRealBasic AS R, LongRealTrans AS RT, RandomIteratedSquaring AS IterSqr, RandomSubtractiveFibo1 AS SubFibo1, RandomSubtractiveFibo2 AS SubFibo2, RandomMultiplicativeFibo AS MulFibo, RandomQuaternaryFibo AS QuaFibo, RandomImprovedMcGill AS McGill, RandomWolframCA AS Wolf, Word, FloatMode; IMPORT RandomRep; <* UNUSED *> CONST Module = "RandomCombinedSlow."; REVEAL T = TPublic BRANDED OBJECT subfibo1: SubFibo1.T; subfibo2: SubFibo2.T; mulfibo : MulFibo.T; quafibo : QuaFibo.T; mcgill : McGill.T; wolf : Wolf.T; OVERRIDES init := Init; generateWord := GenerateWord; generateReal := GenerateReal; END;***************************************************** The random words output by this generator ought to be extremely random since this is a combination of 5 generators, each pretty good by itSELF, and all 5 work according to different principles. Its only disadvantage is it is too slow for applications which do a very small amount of computing per random number. On a 100 MHz Pentium, I generated 10^7 random words in 47 seconds, i.e. 4.7 microseconds per rand, i.e. (assuming roughly 1 instruction per clock) 470 instructions are needed to generate a rand, via RandWord(). Rough number of pentium clocks per call for various routines using the parameters in the CONST declarations (Optimizer on. Asserts on.): SubtractiveFibo2 35 MultiplicativeFibo1 40 QuaternaryFibo 50 SubtractiveFibo1 60 (code looks the fastest, but LONGREALs not Word.T) ImprovedMcGill 100 WolframCA 130 (with wolfnum=5) RandWord 470 (yes, I know the sum isn't 470; timings crude!) Uni01 1000 FasterUni01 80 FasterRandWord 80 Math.sin 130 (for comparison) *********************************************
PROCEDURE** Initializes all random number generators here. Quite slow. If fixed=FALSE (the default) will use the time as seed. If TRUE will use a particular fixed seed. ************************************************************
GenerateWord(SELF: T; ): Word.T = BEGIN RETURN Word.Plus( Word.Plus(Word.Plus(SELF.subfibo2.engine(), SELF.mulfibo.engine()), Word.Plus(SELF.quafibo.engine(), SELF.mcgill.engine())), ORD(SELF.wolf.engine())); END GenerateWord; PROCEDURE GenerateReal(SELF: T; ): R.T = <* FATAL FloatMode.Trap *> VAR x: R.T; BEGIN x := R.Scalb(R.Scalb(FLOAT(GenerateWord(SELF), R.T), 6 - Word.Size) + FLOAT(GenerateWord(SELF), R.T), -Word.Size); <* ASSERT -RT.Half <= x *> <* ASSERT x < 0.52D0 *> IF x < R.Zero THEN x := x + R.One; END; x := x - SELF.subfibo1.engine(); IF x < R.Zero THEN x := x + R.One; END; <* ASSERT x >= R.Zero *> <* ASSERT x < R.One *> RETURN x; END GenerateReal;
Init(SELF: T; fixed: BOOLEAN := FALSE; ): T = BEGIN WITH is = NEW(IterSqr.T).init(fixed) DO SELF.subfibo1 := NEW(SubFibo1.T).init(is); SELF.subfibo2 := NEW(SubFibo2.T).init(is); SELF.mulfibo := NEW(MulFibo.T).init(is); SELF.quafibo := NEW(QuaFibo.T).init(is); SELF.mcgill := NEW(McGill.T).init(is); SELF.wolf := NEW(Wolf.T).init(is); END; (* rev 'em up by 60 calls to Uni01() *) FOR i := 0 TO 60 DO EVAL SELF.generateReal(); END; RETURN SELF; END Init; BEGIN END RandomCombinedSlow.