c7-1 (779515), страница 4
Текст из файла (страница 4)
The resulting in-line generator (we will call it ranqd1) is simplyunsigned long idum;...idum = 1664525L*idum + 1013904223L;This is about as good as any 32-bit linear congruential generator, entirely adequate for manyuses. 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 cangenerate 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 byfloating-point 232 , a dirty trick is to mask in an exponent that makes the value lie between 1 and2, then subtract 1.0. The resulting in-line generator (call it ranqd2) will look something likeSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).unsigned long jran,ia,ic,im;float ran;...jran=(jran*ia+ic) % im;ran=(float) jran / (float) im;2857.1 Uniform DeviatesConstants for Quick and Dirty Random Number Generatorsoverflow atiaic6075106128378752111663787542116636075 13666655 93611979 430128313992531overflow at220221iaic86436 1093121500 1021259200 421182572567354773117128 1277121500 2041312500 74124749256736603714580017500023328024494430809369794929751749227222228223144062928253125967 3041419 6173171 11213212960 1741 273114000 1541 295721870 1291 462131104 625 6571139968 205 29573225366126611861159722924226im139968 3877 29573214326 3613 45289714025 1366 150889230134456 8121259200 7141284115477323129282 1255 617381000 421 17117134456 281 28411233280 9301 49297714025 4096 150889232unsigned long idum,itemp;float rand;#ifdef vaxstatic unsigned long jflone = 0x00004080;static unsigned long jflmsk = 0xffff007f;#elsestatic unsigned long jflone = 0x3f800000;static unsigned long jflmsk = 0x007fffff;#endif...idum = 1664525L*idum + 1013904223L;itemp = jflone | (jflmsk & idum);rand = (*(float *)&itemp)-1.0;The hex constants 3F800000 and 007FFFFF are the appropriate ones for computers usingthe IEEE representation for 32-bit floating-point numbers (e.g., IBM PCs and most UNIXworkstations).
For DEC VAXes, the correct hex constants are, respectively, 00004080 andFFFF007F. Notice that the IEEE mask results in the floating-point number being constructedout of the 23 low-order bits of the integer, which is not ideal. (Your authors have triedvery hard to make almost all of the material in this book machine and compiler independent— indeed, even programming language independent. This subsection is a rare aberration.Forgive us. Once in a great while the temptation to be really dirty is just irresistible.)Relative Timings and RecommendationsTimings are inevitably machine dependent. Nevertheless the following tableis indicative of the relative timings, for typical machines, of the various uniformSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).im286Chapter 7.Random Numbersgenerators discussed in this section, plus ran4 from §7.5. Smaller values in the tableindicate faster generators. The generators ranqd1 and ranqd2 refer to the “quickand dirty” generators immediately above.Relative Execution Timeran0≡ 1.0ran1ran2ran3≈ 1.3≈ 2.0≈ 0.6ranqd1ranqd2ran4≈ 0.10≈ 0.25≈ 4.0On balance, we recommend ran1 for general use. It is portable, based onPark and Miller’s Minimal Standard generator with an additional shuffle, and has noknown (to us) flaws other than period exhaustion.If you are generating more than 100,000,000 random numbers in a singlecalculation (that is, more than about 5% of ran1’s period), we recommend the useof ran2, with its much longer period.Knuth’s subtractive routine ran3 seems to be the timing winner among portableroutines.
Unfortunately the subtractive method is not so well studied, and not astandard. We like to keep ran3 in reserve for a “second opinion,” substituting it whenwe suspect another generator of introducing unwanted correlations into a calculation.The routine ran4 generates extremely good random deviates, and has someother nice properties, but it is slow. See §7.5 for discussion.Finally, the quick and dirty in-line generators ranqd1 and ranqd2 are very fast,but they are somewhat machine dependent, and at best only as good as a 32-bit linearcongruential generator ever is — in our view not good enough in many situations.We would use these only in very special cases, where speed is critical.CITED REFERENCES AND FURTHER READING:Park, S.K., and Miller, K.W.
1988, Communications of the ACM, vol. 31, pp. 1192–1201. [1]Schrage, L. 1979, ACM Transactions on Mathematical Software, vol. 5, pp. 132–138. [2]Bratley, P., Fox, B.L., and Schrage, E.L. 1983, A Guide to Simulation (New York: SpringerVerlag). [3]Knuth, D.E. 1981, Seminumerical Algorithms, 2nd ed., vol.
2 of The Art of Computer Programming(Reading, MA: Addison-Wesley), §§3.2–3.3. [4]Kahaner, D., Moler, C., and Nash, S. 1989, Numerical Methods and Software (Englewood Cliffs,NJ: Prentice Hall), Chapter 10. [5]L’Ecuyer, P. 1988, Communications of the ACM, vol. 31, pp. 742–774. [6]Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for MathematicalComputations (Englewood Cliffs, NJ: Prentice-Hall), Chapter 10.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Generator2877.2 Transformation Method: Exponential and Normal Deviates7.2 Transformation Method: Exponential andNormal DeviatesThe probability distribution p(x) is of course normalized, so thatZ ∞p(x)dx = 1(7.2.2)−∞Now suppose that we generate a uniform deviate x and then take some prescribedfunction of it, y(x). The probability distribution of y, denoted p(y)dy, is determinedby the fundamental transformation law of probabilities, which is simply|p(y)dy| = |p(x)dx| dx p(y) = p(x) dyor(7.2.3)(7.2.4)Exponential DeviatesAs an example, suppose that y(x) ≡ − ln(x), and that p(x) is as given byequation (7.2.1) for a uniform deviate.
Then dx (7.2.5)p(y)dy = dy = e−y dydywhich is distributed exponentially. This exponential distribution occurs frequentlyin real problems, usually as the distribution of waiting times between independentPoisson-random events, for example the radioactive decay of nuclei. You can alsoeasily see (from 7.2.4) that the quantity y/λ has the probability distribution λe−λy .So we have#include <math.h>float expdev(long *idum)Returns an exponentially distributed, positive,ran1(idum) as the source of uniform deviates.{float ran1(long *idum);float dum;dodum=ran1(idum);while (dum == 0.0);return -log(dum);}random deviate of unit mean,usingSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).In the previous section, we learned how to generate random deviates witha uniform probability distribution, so that the probability of generating a numberbetween x and x + dx, denoted p(x)dx, is given byndx 0 < x < 1(7.2.1)p(x)dx =0otherwise.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
