7.1 Uniform Deviates 275 As for references on this subject,the one to turn to first is Knuth [1].Then try [2].Only a few of the standard books on numerical methods [3-4]treat topics relating to random numbers. CITED REFERENCES AND FURTHER READING: Knuth,D.E.1981,Seminumerical Algorithms,2nd ed..vol.2 of The Art of Computer Programming (Reading.MA:Addison-Wesley),Chapter 3,especially $3.5.[1] Bratley.P..Fox.B.L..and Schrage,E.L.1983.A Guide to Simulation (New York:Springer- Verlag).[2] Dahlquist,G.,and Bjorck,A.1974,Numerica/Methods (Englewood Cliffs,NJ:Prentice-Hall). Chapter 11.[3] Forsythe,G.E.,Malcolm,M.A.,and Moler,C.B.1977,Computer Methods for Mathematical Computations(Englewood Cliffs,NJ:Prentice-Hall),Chapter 10.[4] 7.1 Uniform Deviates Uniform deviates are just random numbers that lie within a specified range (typically 0 to 1),with any one number in the range just as likely as any other.They are,in other words,what you probably think "random numbers"are.However, we want to distinguish uniform deviates from other sorts of random numbers,for example numbers drawn from a normal(Gaussian)distribution of specified mean and standard deviation.These other sorts of deviates are almost always generated by performing appropriate operations on one or more uniform deviates,as we will see Programs in subsequent sections.So,a reliable source of random uniform deviates,the subject of this section,is an essential building block for any sort of stochastic modeling or OF SCIENTIFIC Monte Carlo computer work. 61 System-Supplied Random Number Generators Most C implementations have,lurking within,a pair of library routines for initializing,and then generating,"random numbers."In ANSI C,the synopsis is: #include #define RAND_MAX... void srand(unsigned seed); Numerical Recipes 10621 43106 int rand(void); You initialize the random number generator by invoking srand(seed)with (outside some arbitrary seed.Each initializing value will typically result in a different North Software. random sequence,or a least a different starting point in some one enormously long sequence.The same initializing value of seed will always return the same random sequence,however. You obtain successive random numbers in the sequence by successive calls to rand().That function returns an integer that is typically in the range 0 to the largest representable positive value of type int (inclusive).Usually,as in ANSIC, this largest value is available as RAND_MAX,but sometimes you have to figure it out for yourself.If you want a random float value between 0.0 (inclusive)and 1.0 (exclusive),you get it by an expression like
7.1 Uniform Deviates 275 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 machinereadable 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). As for references on this subject, the one to turn to first is Knuth [1]. Then try [2]. Only a few of the standard books on numerical methods [3-4] treat topics relating to random numbers. CITED REFERENCES AND FURTHER READING: Knuth, D.E. 1981, Seminumerical Algorithms, 2nd ed., vol. 2 of The Art of Computer Programming (Reading, MA: Addison-Wesley), Chapter 3, especially §3.5. [1] Bratley, P., Fox, B.L., and Schrage, E.L. 1983, A Guide to Simulation (New York: SpringerVerlag). [2] Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), Chapter 11. [3] Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), Chapter 10. [4] 7.1 Uniform Deviates Uniform deviates are just random numbers that lie within a specified range (typically 0 to 1), with any one number in the range just as likely as any other. They are, in other words, what you probably think “random numbers” are. However, we want to distinguish uniform deviates from other sorts of random numbers, for example numbers drawn from a normal (Gaussian) distribution of specified mean and standard deviation. These other sorts of deviates are almost always generated by performing appropriate operations on one or more uniform deviates, as we will see in subsequent sections. So, a reliable source of random uniform deviates, the subject of this section, is an essential building block for any sort of stochastic modeling or Monte Carlo computer work. System-Supplied Random Number Generators Most C implementations have, lurking within, a pair of library routines for initializing, and then generating, “random numbers.” In ANSI C, the synopsis is: #include #define RAND_MAX ... void srand(unsigned seed); int rand(void); You initialize the random number generator by invoking srand(seed) with some arbitrary seed. Each initializing value will typically result in a different random sequence, or a least a different starting point in some one enormously long sequence. The same initializing value of seed will always return the same random sequence, however. You obtain successive random numbers in the sequence by successive calls to rand(). That function returns an integer that is typically in the range 0 to the largest representable positive value of type int (inclusive). Usually, as in ANSI C, this largest value is available as RAND_MAX, but sometimes you have to figure it out for yourself. If you want a random float value between 0.0 (inclusive) and 1.0 (exclusive), you get it by an expression like
276 Chapter 7.Random Numbers x rand()/(RAND_MAX+1.0); Now our first,and perhaps most important,lesson in this chapter is:be very. very suspicious of a system-supplied rand()that resembles the one just described. If all scientific papers whose results are in doubt because of bad rand()s were to disappear from library shelves,there would be a gap on each shelf about as big as your fist.System-supplied rand()s are almost always linear congruential generators,which generate a sequence of integers I1,I2,13,...,each between 0 and m -1 (e.g.,RAND_MAX)by the recurrence relation 1j+1=ali+c (mod m) (7.1.1) Here m is called the modulus,and a and c are positive integers called the multiplier and the increment respectively.The recurrence (7.1.1)will eventually repeat itself. with a period that is obviously no greater than m.Ifm,a,and c are properly chosen, then the period will be of maximal length,i.e.,of length m.In that case,all possible g integers between 0 and m-1 occur at some point,so any initial"seed"choice of Io is as good as any other:the sequence just takes off from that point. Although this general framework is powerful enough to provide quite decent random numbers,its implementation in many,if not most,ANSI C libraries is quite flawed;quite a number of implementations are in the category "totally botched." Blame should be apportioned about equally between the ANSI C committee and the implementors.The typical problems are these:First,since the ANSI standard specifies that rand()return a value of type int-which is only a two-byte quantity on many machines-RAND_MAX is often not very large.The ANSI C standard requires only that it be at least 32767.This can be disastrous in many circumstances: for a Monte Carlo integration(87.6 and $7.8),you might well want to evaluate 106 N兰 to dir different points,but actually be evaluating the same 32767 points 30 times each,not at all the same thing!You should categorically reject any library random number routine with a two-byte returned value. Second,the ANSI committee's published rationale includes the following mischievous passage:"The committee decided that an implementation should be allowed to provide a rand function which generates the best random sequence 10621 possible in that implementation,and therefore mandated no standard algorithm.It Numerical Recipes 43106 recognized the value,however,of being able to generate the same pseudo-random sequence in different implementations,and so it has published an example.... [emphasis added'The“example'”is (outside Software. unsigned long next=1; North int rand(void)/*NOT RECOMMENDED (see text)*/ next=next*1103515245+12345; return (unsigned int)(next/65536)%32768; void srand(unsigned int seed) next=seed;
276 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 machinereadable 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). x = rand()/(RAND_MAX+1.0); Now our first, and perhaps most important, lesson in this chapter is: be very, very suspicious of a system-supplied rand() that resembles the one just described. If all scientific papers whose results are in doubt because of bad rand()s were to disappear from library shelves, there would be a gap on each shelf about as big as your fist. System-supplied rand()s are almost always linear congruential generators, which generate a sequence of integers I1, I2, I3,..., each between 0 and m − 1 (e.g., RAND_MAX) by the recurrence relation Ij+1 = aIj + c (mod m) (7.1.1) Here m is called the modulus, and a and c are positive integers called the multiplier and the increment respectively. The recurrence (7.1.1) will eventually repeat itself, with a period that is obviously no greater than m. If m, a, and c are properly chosen, then the period will be of maximal length, i.e., of length m. In that case, all possible integers between 0 and m − 1 occur at some point, so any initial “seed” choice of I 0 is as good as any other: the sequence just takes off from that point. Although this general framework is powerful enough to provide quite decent random numbers, its implementation in many, if not most, ANSI C libraries is quite flawed; quite a number of implementations are in the category “totally botched.” Blame should be apportioned about equally between the ANSI C committee and the implementors. The typical problems are these: First, since the ANSI standard specifies that rand() return a value of type int — which is only a two-byte quantity on many machines — RAND_MAX is often not very large. The ANSI C standard requires only that it be at least 32767. This can be disastrous in many circumstances: for a Monte Carlo integration (§7.6 and §7.8), you might well want to evaluate 10 6 different points, but actually be evaluating the same 32767 points 30 times each, not at all the same thing! You should categorically reject any library random number routine with a two-byte returned value. Second, the ANSI committee’s published rationale includes the following mischievous passage: “The committee decided that an implementation should be allowed to provide a rand function which generates the best random sequence possible in that implementation, and therefore mandated no standard algorithm. It recognized the value, however, of being able to generate the same pseudo-random sequence in different implementations, and so it has published an example... . [emphasis added]” The “example” is unsigned long next=1; int rand(void) /* NOT RECOMMENDED (see text) */ { next = next*1103515245 + 12345; return (unsigned int)(next/65536) % 32768; } void srand(unsigned int seed) { next=seed; }
7.1 Uniform Deviates 277 This corresponds to equation (7.1.1)with a 1103515245.c 12345,and m =232(since arithmetic done on unsigned long quantities is guaranteed to return the correct low-order bits).These are not particularly good choices for a and c(the period is only 230),though they are not gross embarrassments by themselves. The real botches occur when implementors,taking the committee's statement above as license,try to "improve"on the published example.For example,one popular 32-bit PC-compatible compiler provides a long generator that uses the above congruence,but swaps the high-order and low-order 16 bits of the returned value. Somebody probably thought that this extra flourish added randomness;in fact it ruins the generator.While these kinds of blunders can,of course,be fixed,there remains a fundamental flaw in simple linear congruential generators,which we now discuss. 8g The linear congruential method has the advantage of being very fast,requiring nted for only a few operations per call,hence its almost universal use.It has the disadvantage g that it is not free of sequential correlation on successive calls.Ifk random numbers at a time are used to plot points in k dimensional space (with each coordinate between 0 and 1),then the points will not tend to "fill up"the k-dimensional space,but rather will lie on (k-1)-dimensional "planes."There will be at most about m/k such planes.If the constants m,a,and c are not very carefully chosen,there will be many fewer than that.If m is as bad as 32768,then the number of planes on which triples of points lie in three-dimensional space will be no greater than about the cube root of 32768.or 32.Even if m is close to the machine's largest representable integer, e.g.,~232,the number of planes on which triples of points lie in three-dimensional space is usually no greater than about the cube root of 232,about 1600.You might Programs well be focusing attention on a physical process that occurs in a small fraction of the total volume,so that the discreteness of the planes can be very pronounced. Even worse,you might be using a generator whose choices of m,a,and c have 是 been botched.One infamous such routine,RANDU,with a =65539 and m =231 兰 to dir was widespread on IBM mainframe computers for many years,and widely copied onto other systems [1].One of us recalls producing a "random"plot with only 11 planes,and being told by his computer center's programming consultant that he had misused the random number generator:"We guarantee that each number is 御你A合过 random individually,but we don't guarantee that more than one of them is random." Figure that out Correlation in k-space is not the only weakness oflinear congruential generators. Numerical Recipes -43108. Such generators often have their low-order(least significant)bits much less random than their high-order bits.If you want to generate a random integer between I and (outside 10,you should always do it using high-order bits,as in Software. j=1+(int)(10.0*rand()/(RAND_MAX+1.0)) Ame ying of and never by anything resembling j=1+(rand()%10); (which uses lower-order bits).Similarly you should never try to take apart a "rand()"number into several supposedly random pieces.Instead use separate calls for every piece
7.1 Uniform Deviates 277 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 machinereadable 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). This corresponds to equation (7.1.1) with a = 1103515245, c = 12345, and m = 232 (since arithmetic done on unsigned long quantities is guaranteed to return the correct low-order bits). These are not particularly good choices for a and c (the period is only 230), though they are not gross embarrassments by themselves. The real botches occur when implementors, taking the committee’s statement above as license, try to “improve” on the published example. For example, one popular 32-bit PC-compatible compiler provides a long generator that uses the above congruence, but swaps the high-order and low-order 16 bits of the returned value. Somebody probably thought that this extra flourish added randomness; in fact it ruins the generator. While these kinds of blunders can, of course, be fixed, there remains a fundamental flaw in simple linear congruential generators, which we now discuss. The linear congruential method has the advantage of being very fast, requiring only a few operations per call, hence its almost universal use. It has the disadvantage that it is not free of sequential correlation on successive calls. If k random numbers at a time are used to plot points in k dimensional space (with each coordinate between 0 and 1), then the points will not tend to “fill up” the k-dimensional space, but rather will lie on (k − 1)-dimensional “planes.” There will be at most about m1/k such planes. If the constants m, a, and c are not very carefully chosen, there will be many fewer than that. If m is as bad as 32768, then the number of planes on which triples of points lie in three-dimensional space will be no greater than about the cube root of 32768, or 32. Even if m is close to the machine’s largest representable integer, e.g., ∼ 232, the number of planes on which triples of points lie in three-dimensional space is usually no greater than about the cube root of 2 32, about 1600. You might well be focusing attention on a physical process that occurs in a small fraction of the total volume, so that the discreteness of the planes can be very pronounced. Even worse, you might be using a generator whose choices of m, a, and c have been botched. One infamous such routine, RANDU, with a = 65539 and m = 2 31, was widespread on IBM mainframe computers for many years, and widely copied onto other systems [1]. One of us recalls producing a “random” plot with only 11 planes, and being told by his computer center’s programming consultant that he had misused the random number generator: “We guarantee that each number is random individually, but we don’t guarantee that more than one of them is random.” Figure that out. Correlation in k-space is not the only weakness of linear congruential generators. Such generators often have their low-order (least significant) bits much less random than their high-order bits. If you want to generate a random integer between 1 and 10, you should always do it using high-order bits, as in j=1+(int) (10.0*rand()/(RAND_MAX+1.0)); and never by anything resembling j=1+(rand() % 10); (which uses lower-order bits). Similarly you should never try to take apart a “rand()” number into several supposedly random pieces. Instead use separate calls for every piece
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 r0. 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 machinereadable 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 multiplicative 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:
7.1 Uniform Deviates 279 #define IA 16807 #define IM 2147483647 #define AM (1.0/IM) #define IQ 127773 #define IR 2836 #define MASK 123459876 float rano(long *idum) "Minimal"random number generator of Park and Miller.Returns a uniform random deviate between 0.0 and 1.0.Set or reset idum to any integer value (except the unlikely value MASK) to initialize the sequence;idum must not be altered between calls for successive deviates in a sequence. long k; float ans; g *idum“=MASK; XORing with MASK allows use of zero and other k=(*idum)/IQ; simple bit patterns for idum. *idum=IA*(*idum-k*IQ)-IR*k: Compute idum=(IA*idum)%IM without over- if (*idum O)*idum +IM; flows by Schrage's method. ans=AM*(*idum) Convert idum to a floating result. *idum MASK; Unmask before return. RECIPES return ans; 婆是 Press. The period of rano is 231-22.1 x 109.A peculiarity of generators of the form (7.1.2)is that the value 0 must never be allowed as the initial seed-it perpetuates itself-and it never occurs for any nonzero initial seed.Experience 9 has shown that users always manage to call random number generators with the seed idum=0.That is why rano performs its exclusive-or with an arbitrary constant both SCIENTIFIC on entry and exit.If you are the first user in history to be proofagainst human error, you can remove the two lines with the A operation. 6 Park and Miller discuss two other multipliers a that can be used with the same m =231-1.These are a =48271 (with g=44488 and r =3399)and a 69621 (with g=30845 and r=23902).These can be substituted in the routine rano if desired;they may be slightly superior to Lewis et al.'s longer-tested values.No values other than these should be used. The routine rano is a Minimal Standard,satisfactory for the majority of Recipes Numerica 10621 applications,but we do not recommend it as the final word on random number 43106 generators.Our reason is precisely the simplicity of the Minimal Standard.It is Recipes not hard to think of situations where successive random numbers might be used in a way that accidentally conflicts with the generation algorithm.For example, since successive numbers differ by a multiple of only 1.6 x 104 out of a modulus of North Software. more than 2 x 109,very small random numbers will tend to be followed by smaller than average values.One time in 105,for example,there will be a value<10-6 returned (as there should be),but this will ahways be followed by a value less than about 0.0168.One can easily think of applications involving rare events where this property would lead to wrong results. There are other,more subtle,serial correlations present in rano.For example, if successive points (Ii,I+1)are binned into a two-dimensional plane for i= 1,2,...,N,then the resulting distribution fails the x2test when N is greater than a few x107,much less than the period m-2.Since low-order serial correlations have historically been such a bugaboo,and since there is a very simple way to remove
7.1 Uniform Deviates 279 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 machinereadable 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). #define IA 16807 #define IM 2147483647 #define AM (1.0/IM) #define IQ 127773 #define IR 2836 #define MASK 123459876 float ran0(long *idum) “Minimal” random number generator of Park and Miller. Returns a uniform random deviate between 0.0 and 1.0. Set or reset idum to any integer value (except the unlikely value MASK) to initialize the sequence; idum must not be altered between calls for successive deviates in a sequence. { long k; float ans; *idum ^= MASK; XORing with MASK allows use of zero and other k=(*idum)/IQ; simple bit patterns for idum. *idum=IA*(*idum-k*IQ)-IR*k; Compute idum=(IA*idum) % IM without overif (*idum < 0) *idum += IM; flows by Schrage’s method. ans=AM*(*idum); Convert idum to a floating result. *idum ^= MASK; Unmask before return. return ans; } The period of ran0 is 231 − 2 ≈ 2.1 × 109. A peculiarity of generators of the form (7.1.2) is that the value 0 must never be allowed as the initial seed — it perpetuates itself — and it never occurs for any nonzero initial seed. Experience has shown that users always manage to call random number generators with the seed idum=0. That is why ran0 performs its exclusive-or with an arbitrary constant both on entry and exit. If you are the first user in history to be proof against human error, you can remove the two lines with the ∧ operation. Park and Miller discuss two other multipliers a that can be used with the same m = 231 − 1. These are a = 48271 (with q = 44488 and r = 3399) and a = 69621 (with q = 30845 and r = 23902). These can be substituted in the routine ran0 if desired; they may be slightly superior to Lewis et al.’s longer-tested values. No values other than these should be used. The routine ran0 is a Minimal Standard, satisfactory for the majority of applications, but we do not recommend it as the final word on random number generators. Our reason is precisely the simplicity of the Minimal Standard. It is not hard to think of situations where successive random numbers might be used in a way that accidentally conflicts with the generation algorithm. For example, since successive numbers differ by a multiple of only 1.6 × 10 4 out of a modulus of more than 2 × 109, very small random numbers will tend to be followed by smaller than average values. One time in 106, for example, there will be a value < 10−6 returned (as there should be), but this will always be followed by a value less than about 0.0168. One can easily think of applications involving rare events where this property would lead to wrong results. There are other, more subtle, serial correlations present in ran0. For example, if successive points (Ii, Ii+1) are binned into a two-dimensional plane for i = 1, 2,...,N, then the resulting distribution fails the χ2 test when N is greater than a few ×107, much less than the period m − 2. Since low-order serial correlations have historically been such a bugaboo, and since there is a very simple way to remove
280 Chapter 7.Random Numbers them,we think that it is prudent to do so. The following routine,ran1,uses the Minimal Standard for its random value, but it shuffles the output to remove low-order serial correlations.A random deviate derived from the jth value in the sequence,Ij,is output not on the jth call,but rather on a randomized later call,j+32 on average.The shuffling algorithm is due to Bays and Durham as described in Knuth [4],and is illustrated in Figure 7.1.1. #define IA 16807 #def1neIM2147483647 #define AM (1.0/IM) #define IQ 127773 #define IR 2836 #define NTAB 32 granted for 19881992 #define NDIV (1+(IM-1)/NTAB) #define EPS 1.2e-7 1.800 #define RNMX (1.0-EPS) float rani(long *idum) "Minimal"random number generator of Park and Miller with Bays-Durham shuffle and added NUMERICAL RECIPES safeguards.Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint values).Call with idum a negative integer to initialize;thereafter,do not alter idum between successive deviates in a sequence.RNMX should approximate the largest floating value that is to make less than 1. f int j; ART long k; static long iy=0; static long iv[NTAB]; 9 Programs float temp; if (*idum =0;j--){ Load the shuffle table (after 8 warm-ups) to dir k=(*idum)/IQ; *idum=IA*(*idum-k*IQ)-IR*k; COMPUTING if (*idum 0)*idum +IM; if (jRNMX)return RNMX;Because users don't expect endpoint values. North Software. else return temp; The routine ran1 passes those statistical tests that rano is known to fail.In fact,we do not know of any statistical test that rani fails to pass,except when the number of calls starts to become on the order of the period m,say >10m/20. For situations when even longer random sequences are needed,L'Ecuyer [6]has given a good way of combining two different sequences with different periods so as to obtain a new sequence whose period is the least common multiple of the two periods.The basic idea is simply to add the two sequences,modulo the modulus of
280 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 machinereadable 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). them, we think that it is prudent to do so. The following routine, ran1, uses the Minimal Standard for its random value, but it shuffles the output to remove low-order serial correlations. A random deviate derived from the jth value in the sequence, I j , is output not on the jth call, but rather on a randomized later call, j + 32 on average. The shuffling algorithm is due to Bays and Durham as described in Knuth [4], and is illustrated in Figure 7.1.1. #define IA 16807 #define IM 2147483647 #define AM (1.0/IM) #define IQ 127773 #define IR 2836 #define NTAB 32 #define NDIV (1+(IM-1)/NTAB) #define EPS 1.2e-7 #define RNMX (1.0-EPS) float ran1(long *idum) “Minimal” random number generator of Park and Miller with Bays-Durham shuffle and added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint values). Call with idum a negative integer to initialize; thereafter, do not alter idum between successive deviates in a sequence. RNMX should approximate the largest floating value that is less than 1. { int j; long k; static long iy=0; static long iv[NTAB]; float temp; if (*idum =0;j--) { Load the shuffle table (after 8 warm-ups). k=(*idum)/IQ; *idum=IA*(*idum-k*IQ)-IR*k; if (*idum RNMX) return RNMX; Because users don’t expect endpoint values. else return temp; } The routine ran1 passes those statistical tests that ran0 is known to fail. In fact, we do not know of any statistical test that ran1 fails to pass, except when the number of calls starts to become on the order of the period m, say > 10 8 ≈ m/20. For situations when even longer random sequences are needed, L’Ecuyer [6] has given a good way of combining two different sequences with different periods so as to obtain a new sequence whose period is the least common multiple of the two periods. The basic idea is simply to add the two sequences, modulo the modulus of
7.1 Uniform Deviates 281 iy 1 ivo OUTPUT RAN 3 ② .com or call (including this one) 19881992 11-800-872 g Cambridge iv3l to any server computer, -7423(North America users to make one paper from NUMERICAL RECIPES IN C: e University Press. THE ART Figure 7.1.1.Shuffling procedure used in ran1 to break up sequential correlations in the Minimal Standard generator.Circled numbers indicate the sequence of events:On each call,the random number 9 Programs in iy is used to choose a random element in the array iv.That element becomes the output random ictly proh number,and also is the next iy.Its spot in iv is refilled from the Minimal Standard routine. either of them (call it m).A trick to avoid an intermediate value that overflows the to dir integer wordsize is to subtract rather than add,and then add back the constant m-1 if the result is 0,so as to wrap around into the desired interval 0,...,m-1. OF SCIENTIFIC COMPUTING(ISBN Notice that it is not necessary that this wrapped subtraction be able to reach 1988-19920 all values 0,...,m-1 from every value of the first sequence.Consider the absurd extreme case where the value subtracted was only between I and 10:The resulting 10-521 sequence would still be no less random than the first sequence by itself.As a practical matter it is only necessary that the second sequence have a range covering substantially all of the range of the first.L'Ecuyer recommends the use of the two Numerical Recipes 43198-5 generators m1=2147483563(with a1=40014,q1=53668,r1=12211)and m2=2147483399(with a2=40692,92=52774,r2=3791).Both moduli (outside are slightly less than231.The periods m1-1=2×3×7×631×81031and North Software. m2-1=2×19×31×1019×1789 share only the factor2,so the period of the combined generator is2.3x 1018.For present computers,period exhaustion is a practical impossibility. visit website Combining the two generators breaks up serial correlations to a considerable machine extent.We nevertheless recommend the additional shuffle that is implemented in the following routine,ran2.We think that,within the limits of its floating-point precision,ran2 provides perfect random numbers,a practical definition of"perfect" is that we will pay $1000 to the first reader who convinces us otherwise (by finding a statistical test that ran2 fails in a nontrivial way,excluding the ordinary limitations of a machine's floating-point representation)
7.1 Uniform Deviates 281 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 machinereadable 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). OUTPUT RAN 1 3 2 iy iv0 iv31 Figure 7.1.1. Shuffling procedure used in ran1 to break up sequential correlations in the Minimal Standard generator. Circled numbers indicate the sequence of events: On each call, the random number in iy is used to choose a random element in the array iv. That element becomes the output random number, and also is the next iy. Its spot in iv is refilled from the Minimal Standard routine. either of them (call it m). A trick to avoid an intermediate value that overflows the integer wordsize is to subtract rather than add, and then add back the constant m − 1 if the result is ≤ 0, so as to wrap around into the desired interval 0,...,m − 1. Notice that it is not necessary that this wrapped subtraction be able to reach all values 0,...,m − 1 from every value of the first sequence. Consider the absurd extreme case where the value subtracted was only between 1 and 10: The resulting sequence would still be no less random than the first sequence by itself. As a practical matter it is only necessary that the second sequence have a range covering substantially all of the range of the first. L’Ecuyer recommends the use of the two generators m1 = 2147483563 (with a1 = 40014, q1 = 53668, r1 = 12211) and m2 = 2147483399 (with a2 = 40692, q2 = 52774, r2 = 3791). Both moduli are slightly less than 231. The periods m1 − 1=2 × 3 × 7 × 631 × 81031 and m2 − 1=2 × 19 × 31 × 1019 × 1789 share only the factor 2, so the period of the combined generator is ≈ 2.3 × 1018. For present computers, period exhaustion is a practical impossibility. Combining the two generators breaks up serial correlations to a considerable extent. We nevertheless recommend the additional shuffle that is implemented in the following routine, ran2. We think that, within the limits of its floating-point precision, ran2 provides perfect random numbers; a practical definition of “perfect” is that we will pay $1000 to the first reader who convinces us otherwise (by finding a statistical test that ran2 fails in a nontrivial way, excluding the ordinary limitations of a machine’s floating-point representation)
282 Chapter 7.Random Numbers #define IM12147483563 #def1neIN22147483399 #define AM (1.0/IM1) #define IMM1 (IM1-1) #define IA1 40014 #define IA2 40692 #define IQ1 53668 #define IQ2 52774 #define IR1 12211 #define IR2 3791 #define NTAB 32 #define NDIV (1+IMM1/NTAB) http://www.nr. #define EPS 1.2e-7 #define RNMX (1.0-EPS) 83 float ran2(long *idum) granted for 19881992 Long period (>2x 1018)random number generator of L'Ecuyer with Bays-Durham shuffle 鱼 and added safeguards.Returns a uniform random deviate between 0.0 and 1.0 (exclusive of 1800 the endpoint values).Call with idum a negative integer to initialize:thereafter,do not alter idum between successive deviates in a sequence.RNMX should approximate the largest floating value that is less than 1. from NUMERICAL RECIPESI int ji long ki static long idum2=123456789; (Nort static long iy=0; static long iv[NTAB]; America server computer, to make e University Press. THE float temp; one paper ART if (*idum =0;j--){ Load the shuffle table (after 8 warm-ups). SCIENTIFIC k=(*idum)/IQ1; *idum=IA1*(*idum-k*IQ1)-k*IR1; to dir if (*idum 0)*idum +IM1; if (j NTAB)iv[j]tidum; 188 iy=iv[o]; 1920 COMPUTING(ISBN k=(*idum)/IQ1; Start here when not initializing. *idum=IA1*(*idum-k*IQ1)-k*IR1; Compute idum=(IA1*idum)%IM1 without 0621 if (*idum 0)*idum +IM1; overflows by Schrage's method. k=idum2/IQ2; idum2=IA2*(idum2-k*IQ2)-k*IR2; Compute idum2=(IA2*idum)%IM2 likewise. Numerical Recipes -43108 if (idum2 0)idum2 +IM2; i=iy/NDIV; Will be in the range 0..NTAB-1. iy=iv[j]-idum2; Here idum is shuffled,idum and idum2 are (outside iv[j]*idum; combined to generate output. if (iy 1)iy +IMM1; North Software. if ((temp=AM*iy)>RNMX)return RNMX; Because users don't expect endpoint values. else return temp; Ame visit website machine L'Ecuyer [6]lists additional short generators that can be combined into longer ones,including generators that can be implemented in 16-bit integer arithmetic. Finally,we give you Knuth's suggestion [4]for a portable routine,which we have translated to the present conventions as ran3.This is not based on the linear congruential method at all,but rather on a subtractive method (see also [5)).One might hope that its weaknesses,if any,are therefore of a highly different character
282 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 machinereadable 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). #define IM1 2147483563 #define IM2 2147483399 #define AM (1.0/IM1) #define IMM1 (IM1-1) #define IA1 40014 #define IA2 40692 #define IQ1 53668 #define IQ2 52774 #define IR1 12211 #define IR2 3791 #define NTAB 32 #define NDIV (1+IMM1/NTAB) #define EPS 1.2e-7 #define RNMX (1.0-EPS) float ran2(long *idum) Long period (> 2 × 1018) random number generator of L’Ecuyer with Bays-Durham shuffle and added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint values). Call with idum a negative integer to initialize; thereafter, do not alter idum between successive deviates in a sequence. RNMX should approximate the largest floating value that is less than 1. { int j; long k; static long idum2=123456789; static long iy=0; static long iv[NTAB]; float temp; if (*idum =0;j--) { Load the shuffle table (after 8 warm-ups). k=(*idum)/IQ1; *idum=IA1*(*idum-k*IQ1)-k*IR1; if (*idum RNMX) return RNMX; Because users don’t expect endpoint values. else return temp; } L’Ecuyer [6] lists additional short generators that can be combined into longer ones, including generators that can be implemented in 16-bit integer arithmetic. Finally, we give you Knuth’s suggestion [4] for a portable routine, which we have translated to the present conventions as ran3. This is not based on the linear congruential method at all, but rather on a subtractive method (see also [5]). One might hope that its weaknesses, if any, are therefore of a highly different character
7.1 Uniform Deviates 283 from the weaknesses,if any,of ran1 above.If you ever suspect trouble with one routine,it is a good idea to try the other in the same application.ran3 has one nice feature:if your machine is poor on integer arithmetic (i.e.,is limited to 16-bit integers),you can declare mj,mk,and ma[]as float,define mbig and mseed as 4000000 and 1618033,respectively,and the routine will be rendered entirely floating-point. #include Change to math.h in K&R C. #define MBIG 1000000000 #define MSEED 161803398 #define MZ 0 #define FAC (1.0/MBIG) According to Knuth,any large MBIG,and any smaller (but still large)MSEED can be substituted 83g (includi for the above values. 18881892 float ran3(long *idum) 11.800 Returns a uniform random deviate between 0.0 and 1.0.Set idum to any negative value to initialize or reinitialize the sequence. f from NUMERICAL RECIPESI static int inext,inextp; 6 static long ma[56]; The value 56 (range ma[1..55])is special and static int iff-0; should not be modified;see Knuth. long mj,mk; int i,ii,k; Press. if (*idum o ll iff ==0){ Initialization ART 1ff=1; mj=labs(MSEED-labs(*idum)); Initialize ma[55]using the seed idum and the 9 mj %MBIG: large number MSEED. Programs 豆 ma[55]-mj; mk=1; SCIENTIFIC for(1=1;1<=54;1++)[ Now initialize the rest of the table, 11=(21*1)%55; in a slightly random order, ma[ii]=mk: with numbers that are not especially random. mk=mj-mk; if (mk MZ)mk +MBIG; mj=ma [ii]; 1920 COMPUTING(ISBN for(k=1;k<=4;k++) We randomize them by "warming up the gener- for(1=1;1<=55;1+)[ ator." ma[1]-=ma[1+(1+30)%55]; if (ma[i]MZ)ma[i]+MBIG; inext=0; Prepare indices for our first generated number Numerical Recipes -43108 inextp=31; The constant 31 is special;see Knuth. *idum=1; (outside Here is where we start,except on initialization. Software. if (++inext =56)inext=1; Increment inext and inextp,wrapping around North if (++inextp ==56)inextp=1; 56to1. mj=ma[inext]-ma[inextp]; Generate a new random number subtractively. Ame if (mj MZ)mj +MBIG; Be sure that it is in range. ma[inext]=mj; Store it, return mj*FAC; and output the derived uniform deviate. Quick and Dirty Generators One sometimes would like a"quick and dirty"generator to embed in a program,perhaps taking only one or two lines of code,just to somewhat randomize things.One might wish to
7.1 Uniform Deviates 283 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 machinereadable 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). from the weaknesses, if any, of ran1 above. If you ever suspect trouble with one routine, it is a good idea to try the other in the same application. ran3 has one nice feature: if your machine is poor on integer arithmetic (i.e., is limited to 16-bit integers), you can declare mj, mk, and ma[] as float, define mbig and mseed as 4000000 and 1618033, respectively, and the routine will be rendered entirely floating-point. #include Change to math.h in K&R C. #define MBIG 1000000000 #define MSEED 161803398 #define MZ 0 #define FAC (1.0/MBIG) According to Knuth, any large MBIG, and any smaller (but still large) MSEED can be substituted for the above values. float ran3(long *idum) Returns a uniform random deviate between 0.0 and 1.0. Set idum to any negative value to initialize or reinitialize the sequence. { static int inext,inextp; static long ma[56]; The value 56 (range ma[1..55]) is special and static int iff=0; should not be modified; see Knuth. long mj,mk; int i,ii,k; if (*idum < 0 || iff == 0) { Initialization. iff=1; mj=labs(MSEED-labs(*idum)); Initialize ma[55] using the seed idum and the mj %= MBIG; large number MSEED. ma[55]=mj; mk=1; for (i=1;i<=54;i++) { Now initialize the rest of the table, ii=(21*i) % 55; in a slightly random order, ma[ii]=mk; with numbers that are not especially random. mk=mj-mk; if (mk < MZ) mk += MBIG; mj=ma[ii]; } for (k=1;k<=4;k++) We randomize them by “warming upthe generfor (i=1;i<=55;i++) { ator.” ma[i] -= ma[1+(i+30) % 55]; if (ma[i] < MZ) ma[i] += MBIG; } inext=0; Prepare indices for our first generated number. inextp=31; The constant 31 is special; see Knuth. *idum=1; } Here is where we start, except on initialization. if (++inext == 56) inext=1; Increment inext and inextp, wrapping around if (++inextp == 56) inextp=1; 56 to 1. mj=ma[inext]-ma[inextp]; Generate a new random number subtractively. if (mj < MZ) mj += MBIG; Be sure that it is in range. ma[inext]=mj; Store it, return mj*FAC; and output the derived uniform deviate. } Quick and Dirty Generators One sometimes would like a “quick and dirty” generator to embed in a program, perhaps taking only one or two lines of code, just to somewhat randomize things. One might wish to
284 Chapter 7.Random Numbers process data from an experiment not always in exactly the same order,for example,so that the first output is more"typical"than might otherwise be the case. For this kind of application,all we really need is a list of"good"choices for m,a,and c in equation(7.1.1).If we don't need a period longer than 10 to 106,say,we can keep the value of(m-1)a+e small enough to avoid overflows that would otherwise mandate the extra complexity of Schrage's method (above).We can thus easily embed in our programs unsigned long jran,ia,ic,im; float ran; iran=(iraniatic)%im; ran=(float)jran (float)im; whenever we want a quick and dirty uniform deviate,or jran=(jran*iatic)%im; NUMERICAL j=jlo+((jhi-jlo+1)*jran)/im; whenever we want an integer between ilo and jhi,inclusive.(In both cases jran was once initialized to any seed value between 0 and im-1.) RECIPES I Be sure to remember,however,that when im is small,the kth root of it,which is the 令 number of planes in k-space,is even smaller!So a quick and dirty generator should never be used to select points in k-space with k>1. With these caveats,some"good"choices for the constants are given in the accompanying table.These constants(i)give a period of maximal length im,and,more important,(ii)pass 2 Knuth's "spectral test"for dimensions 2,3,4,5,and 6.The increment ic is a prime,close to the value (-v3)im;actually almost any value of ic that is relatively prime to im will do just as well,but there is some "lore"favoring this choice (see [41,p.84). IENTIFIC An Even Quicker Generator 61 In C,if you multiply two unsigned long int integers on a machine with a 32-bit long integer representation,the value returned is the low-order 32 bits of the true 64-bit product.If we now choose m =232,the "mod"in equation (7.1.1)is free,and we have simply 1j+1=alj+c (7.1.6) 3、 10621 Knuth suggests a =1664525 as a suitable multiplier for this value of m.H.W.Lewis Numerica has conducted extensive tests of this value of a with c=1013904223,which is a prime close 43126 to (v5-2)m.The resulting in-line generator (we will call it ranqd1)is simply E喜 unsigned long idum; 1dum=1664525L*1dum+1013904223L; North This is about as good as any 32-bit linear congruential generator,entirely adequate for many uses.And,with only a single multiply and add,it is very fast. To check whether your machine has the desired integer properties,see if you can generate the following sequence of 32-bit values (given here in hex):00000000,3C6EF35F, 47502932,D1CCF6E9,AAF95334,6252E503,9F2EC686,57FE6C2D,A3D95FA8,81FD- BEE7,94FOAF1A,CBF633B1. If you need floating-point values instead of 32-bit integers,and want to avoid a divide by floating-point 232,a dirty trick is to mask in an exponent that makes the value lie between I and 2,then subtract 1.0.The resulting in-line generator(call it ranqd2)will look something like
284 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 machinereadable 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). process data from an experiment not always in exactly the same order, for example, so that the first output is more “typical” than might otherwise be the case. For this kind of application, all we really need is a list of “good” choices for m, a, and c in equation (7.1.1). If we don’t need a period longer than 104 to 106, say, we can keep the value of (m − 1)a + c small enough to avoid overflows that would otherwise mandate the extra complexity of Schrage’s method (above). We can thus easily embed in our programs unsigned long jran,ia,ic,im; float ran; ... jran=(jran*ia+ic) % im; ran=(float) jran / (float) im; whenever we want a quick and dirty uniform deviate, or jran=(jran*ia+ic) % im; j=jlo+((jhi-jlo+1)*jran)/im; whenever we want an integer between jlo and jhi, inclusive. (In both cases jran was once initialized to any seed value between 0 and im-1.) Be sure to remember, however, that when im is small, the kth root of it, which is the number of planes in k-space, is even smaller! So a quick and dirty generator should never be used to select points in k-space with k > 1. With these caveats, some “good” choices for the constants are given in the accompanying table. These constants (i) give a period of maximal length im, and, more important, (ii) pass Knuth’s “spectral test” for dimensions 2, 3, 4, 5, and 6. The increment ic is a prime, close to the value ( 1 2 − 1 6 √3)im; actually almost any value of ic that is relatively prime to im will do just as well, but there is some “lore” favoring this choice (see [4], p. 84). An Even Quicker Generator In C, if you multiply two unsigned long int integers on a machine with a 32-bit long integer representation, the value returned is the low-order 32 bits of the true 64-bit product. If we now choose m = 232, the “mod” in equation (7.1.1) is free, and we have simply Ij+1 = aIj + c (7.1.6) Knuth suggests a = 1664525 as a suitable multiplier for this value of m. H.W. Lewis has conducted extensive tests of this value of a with c = 1013904223, which is a prime close to ( √5 − 2)m. The resulting in-line generator (we will call it ranqd1) is simply unsigned long idum; ... idum = 1664525L*idum + 1013904223L; This is about as good as any 32-bit linear congruential generator, entirely adequate for many uses. And, with only a single multiply and add, it is very fast. To check whether your machine has the desired integer properties, see if you can generate the following sequence of 32-bit values (given here in hex): 00000000, 3C6EF35F, 47502932, D1CCF6E9, AAF95334, 6252E503, 9F2EC686, 57FE6C2D, A3D95FA8, 81FDBEE7, 94F0AF1A, CBF633B1. If you need floating-point values instead of 32-bit integers, and want to avoid a divide by floating-point 232, a dirty trick is to mask in an exponent that makes the value lie between 1 and 2, then subtract 1.0. The resulting in-line generator (call it ranqd2) will look something like