278 Chapter 7.Random Numbers Portable Random Number Generators Park and Miller [1]have surveyed a large number of random number generators that have been used over the last 30 years or more.Along with a good theoretical review,they present an anecdotal sampling of a number of inadequate generators that have come into widespread use.The historical record is nothing if not appalling. There is good evidence,both theoretical and empirical,that the simple multi- plicative congruential algorithm Ij+1=alj (mod m) (7.1.2) can be as good as any of the more general linear congruential generators that have c0(equation 7.1.1)-if the multiplier a and modulus m are chosen exquisitely carefully.Park and Miller propose a "Minimal Standard"generator based on the choices RECIPES a=75=16807m=231-1=2147483647 (7.1.3) 9 First proposed by Lewis,Goodman,and Miller in 1969,this generator has in subsequent years passed all new theoretical tests,and(perhaps more importantly) has accumulated a large amount of successful use.Park and Miller do not claim that the generator is"perfect"(we will see below that it is not),but only that it is a good minimal standard against which other generators should be judged. It is not possible to implement equations (7.1.2)and (7.1.3)directly in a high-level language,since the product of a and m-1 exceeds the maximum value IENTIFIC for a 32-bit integer.Assembly language implementation using a 64-bit product 6 register is straightforward,but not portable from machine to machine. A trick due to Schrage [2.3]for multiplying two 32-bit integers modulo a 32-bit constant, without using any intermediates larger than 32 bits (including a sign bit)is therefore extremely interesting:It allows the Minimal Standard generator to be implemented in essentially any programming language on essentially any machine. Schrage's algorithm is based on an approximate factorization of m, Numerica 10.621 43126 m=aq+r,i.e.,q=[m/al,r =m mod a (7.1.4) with square brackets denoting integer part.If r is small,specifically r<q,and 腿 0<z<m-1,it can be shown that both a(z mod g)and r[z/g]lie in the range North 0,...,m-1,and that ∫a(z mod q)-rlz/g ifit is >0. az mod m= a(z mod q)-r[z/g]m otherwise (7.1.5) The application of Schrage's algorithm to the constants(7.1.3)uses the values q=127773andr=2836. Here is an implementation of the Minimal Standard generator:278 Chapter 7. Random Numbers Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Portable Random Number Generators Park and Miller [1] have surveyed a large number of random number generators that have been used over the last 30 years or more. Along with a good theoretical review, they present an anecdotal sampling of a number of inadequate generators that have come into widespread use. The historical record is nothing if not appalling. There is good evidence, both theoretical and empirical, that the simple multi￾plicative congruential algorithm Ij+1 = aIj (mod m) (7.1.2) can be as good as any of the more general linear congruential generators that have c = 0 (equation 7.1.1) — if the multiplier a and modulus m are chosen exquisitely carefully. Park and Miller propose a “Minimal Standard” generator based on the choices a = 75 = 16807 m = 231 − 1 = 2147483647 (7.1.3) First proposed by Lewis, Goodman, and Miller in 1969, this generator has in subsequent years passed all new theoretical tests, and (perhaps more importantly) has accumulated a large amount of successful use. Park and Miller do not claim that the generator is “perfect” (we will see below that it is not), but only that it is a good minimal standard against which other generators should be judged. It is not possible to implement equations (7.1.2) and (7.1.3) directly in a high-level language, since the product of a and m − 1 exceeds the maximum value for a 32-bit integer. Assembly language implementation using a 64-bit product register is straightforward, but not portable from machine to machine. A trick due to Schrage [2,3] for multiplying two 32-bit integers modulo a 32-bit constant, without using any intermediates larger than 32 bits (including a sign bit) is therefore extremely interesting: It allows the Minimal Standard generator to be implemented in essentially any programming language on essentially any machine. Schrage’s algorithm is based on an approximate factorization of m, m = aq + r, i.e., q = [m/a], r = m mod a (7.1.4) with square brackets denoting integer part. If r is small, specifically r<q, and 0 <z<m − 1, it can be shown that both a(z mod q) and r[z/q] lie in the range 0,...,m − 1, and that az mod m =
a(z mod q) − r[z/q] if it is ≥ 0, a(z mod q) − r[z/q] + m otherwise (7.1.5) The application of Schrage’s algorithm to the constants (7.1.3) uses the values q = 127773 and r = 2836. Here is an implementation of the Minimal Standard generator: