296 Chapter 7.Random Numbers if (n !nold){ If n has changed,then compute useful quanti- en=n; ties. oldg=gammln(en+1.0); nold=n; if (p !pold){ If p has changed,then compute useful quanti- pc=1.0-p; ties. plog=log(p); pclog=log(pc); pold=p; sq=sqrt(2.0*am*pc); The following code should by now seem familiar: do rejection method with a Lorentzian compar- do ison function. angle=PI*ran1(idum); y=tan(angle); em=sq*y+am; while (em 0.0 II em >(en+1.0)); Reject. em=floor(em); Trick for integer-valued distribution. t=1.2*sq*(1.0+y*y)*exp(oldg-gammln(em+1.0) -gammln(en-em+1.0)+em*plog+(en-em)*pclog); while (ran1(idum)>t); Reject.This happens about 1.5 times per devi- ate,on average if (p !pp)bnl=n-bnl; Remember to undo the symmetry transforma- A 令 return bnl; tion. Press. See Devroye [2]and Bratley [3]for many additional algorithms. CITED REFERENCES AND FURTHER READING: Knuth,D.E.1981.Seminumerical Algorithms,2nded.,vol.2 of The Art of Computer Programming IENTIFIC (Reading,MA:Addison-Wesley),pp.120ff.[1] 61 Devroye,L.1986,Non-Uniform Random Variate Generation (New York:Springer-Verlag),SX.4. 2 Bratley,P.,Fox,B.L.,and Schrage,E.L.1983,A Guide to Simulation (New York:Springer- Verlag）.[3l. COMPUTING (ISBN 188810920 Numerical Recipes 10.621 43108 7.4 Generation of Random Bits (outside The C language gives you useful access to some machine-level bitwise operations North Software. such as<<(left shift).This section will show you how to put such abilities to good use. The problem is how to generate single random bits,with 0 and 1 equally probable.Of course you can just generate uniform random deviates between zero and one and use their high-order bit(i.e.,test if they are greater than or less than 0.5).However this takes a lot of arithmetic;there are special-purpose applications, such as real-time signal processing,where you want to generate bits very much faster than that. One method for generating random bits,with two variant implementations,is based on "primitive polynomials modulo 2."The theory of these polynomials is beyond our scope (although 87.7 and $20.3 will give you small tastes of it).Here,296 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). if (n != nold) { If n has changed, then compute useful quanti￾en=n; ties. oldg=gammln(en+1.0); nold=n; } if (p != pold) { If p has changed, then compute useful quanti￾pc=1.0-p; ties. plog=log(p); pclog=log(pc); pold=p; } sq=sqrt(2.0*am*pc); The following code should by now seem familiar: rejection method with a Lorentzian compar￾ison function. do { do { angle=PI*ran1(idum); y=tan(angle); em=sq*y+am; } while (em < 0.0 || em >= (en+1.0)); Reject. em=floor(em); Trick for integer-valued distribution. t=1.2*sq*(1.0+y*y)*exp(oldg-gammln(em+1.0) -gammln(en-em+1.0)+em*plog+(en-em)*pclog); } while (ran1(idum) > t); Reject. This happens about 1.5 times per devi￾bnl=em; ate, on average. } if (p != pp) bnl=n-bnl; Remember to undo the symmetry transforma￾return bnl; tion. } See Devroye [2] and Bratley [3] for many additional algorithms. 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), pp. 120ff. [1] Devroye, L. 1986, Non-Uniform Random Variate Generation (New York: Springer-Verlag), §X.4. [2] Bratley, P., Fox, B.L., and Schrage, E.L. 1983, A Guide to Simulation (New York: Springer￾Verlag). [3]. 7.4 Generation of Random Bits The C language gives you useful access to some machine-level bitwise operations such as << (left shift). This section will show you how to put such abilities to good use. The problem is how to generate single random bits, with 0 and 1 equally probable. Of course you can just generate uniform random deviates between zero and one and use their high-order bit (i.e., test if they are greater than or less than 0.5). However this takes a lot of arithmetic; there are special-purpose applications, such as real-time signal processing, where you want to generate bits very much faster than that. One method for generating random bits, with two variant implementations, is based on “primitive polynomials modulo 2.” The theory of these polynomials is beyond our scope (although §7.7 and §20.3 will give you small tastes of it). Here