Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 73
Текст из файла (страница 73)
For greater portability, you can insteadconstruct a floating number by making the (signed) 32-bit integer nonnegative (typically, youadd exactly 231 if it is negative) and then multiplying it by a floating constant (typically 2.−31 ).An interesting, and sometimes useful, feature of the routine ran4, below, is that it allowsrandom access to the nth random value in a sequence, without the necessity of first generatingvalues 1 · · · n − 1.
This property is shared by any random number generator based on hashing(the technique of mapping data keys, which may be highly clustered in value, approximatelyuniformly into a storage address space) [5,6] . One might have a simulation problem in whichsome certain rare situation becomes recognizable by its consequences only considerably afterit has occurred. One may wish to restart the simulation back at that occurrence, using identicalrandom values but, say, varying some other control parameters.
The relevant question mightthen be something like “what random numbers were used in cycle number 337098901?” Itmight already be cycle number 395100273 before the question comes up. Random generatorsbased on recursion, rather than hashing, cannot easily answer such a question.float ran4(long *idum)Returns a uniform random deviate in the range 0.0 to 1.0, generated by pseudo-DES (DESlike) hashing of the 64-bit word (idums,idum), where idums was set by a previous call withnegative idum. Also increments idum. Routine can be used to generate a random sequenceby successive calls, leaving idum unaltered between calls; or it can randomly access the nthdeviate in a sequence by calling with idum = n.
Different sequences are initialized by calls withdiffering negative values of idum.{void psdes(unsigned long *lword, unsigned long *irword);unsigned long irword,itemp,lword;static long idums = 0;The hexadecimal constants jflone and jflmsk below are used to produce a floating numberbetween 1.
and 2. by bitwise masking. They are machine-dependent. See text.#if defined(vax) || defined(_vax_) || defined(__vax__) || defined(VAX)static unsigned long jflone = 0x00004080;static unsigned long jflmsk = 0xffff007f;#elsestatic unsigned long jflone = 0x3f800000;static unsigned long jflmsk = 0x007fffff;#endifif (*idum < 0) {idums = -(*idum);*idum=1;}irword=(*idum);lword=idums;Reset idums and prepare to return the firstdeviate in its sequence.304Chapter 7.psdes(&lword,&irword);itemp=jflone | (jflmsk & irword);++(*idum);return (*(float *)&itemp)-1.0;Random Numbers“Pseudo-DES” encode the words.Mask to a floating number between 1 and2.Subtraction moves range to 0.
to 1.}The accompanying table gives data for verifying that ran4 and psdes work correctlyon your machine. We do not advise the use of ran4 unless you are able to reproduce thehex values shown. Typically, ran4 is about 4 times slower than ran0 (§7.1), or about 3times slower than ran1.Values for Verifying the Implementation of psdesidumbefore psdes callafter psdes call (hex)ran4(idum)lwordirwordlwordirwordVAXPC–111604D1DCE509C0C230.2758980.21912099199D97F8571A66CB41A0.2082040.849246–999917822309D643009840.0343070.375290999999D7F376F059BA89EB0.8386760.457334Successive calls to psdes with arguments −1, 99, −99, and 1, should produce exactly thelword and irword values shown. Masking conversion to a returned floating random valueis allowed to be machine dependent; values for VAX and PC are shown.CITED REFERENCES AND FURTHER READING:Data Encryption Standard, 1977 January 15, Federal Information Processing Standards Publication, number 46 (Washington: U.S.
Department of Commerce, National Bureau of Standards). [1]Guidelines for Implementing and Using the NBS Data Encryption Standard, 1981 April 1, FederalInformation Processing Standards Publication, number 74 (Washington: U.S. Departmentof Commerce, National Bureau of Standards). [2]Validating the Correctness of Hardware Implementations of the NBS Data Encryption Standard,1980, NBS Special Publication 500–20 (Washington: U.S. Department of Commerce, National Bureau of Standards). [3]Meyer, C.H. and Matyas, S.M. 1982, Cryptography: A New Dimension in Computer Data Security(New York: Wiley).
[4]Knuth, D.E. 1973, Sorting and Searching, vol. 3 of The Art of Computer Programming (Reading,MA: Addison-Wesley), Chapter 6. [5]Vitter, J.S., and Chen, W-C. 1987, Design and Analysis of Coalesced Hashing (New York:Oxford University Press). [6]7.6 Simple Monte Carlo IntegrationInspirations for numerical methods can spring from unlikely sources. “Splines”first were flexible strips of wood used by draftsmen. “Simulated annealing” (weshall see in §10.9) is rooted in a thermodynamic analogy. And who does not feel atleast a faint echo of glamor in the name “Monte Carlo method”?7.6 Simple Monte Carlo Integration305Suppose that we pick N random points, uniformly distributed in a multidimensional volume V . Call them x1 , . .
. , xN . Then the basic theorem of Monte Carlointegration estimates the integral of a function f over the multidimensional volume,*&f 2 − f2(7.6.1)f dV ≈ V f ± VNHere the angle brackets denote taking the arithmetic mean over the N sample points,f ≡N1 f(xi )N i=11N21 2f2 ≡f (xi )N i=1(7.6.2)The “plus-or-minus” term in (7.6.1) is a one standard deviation error estimate forthe integral, not a rigorous bound; further, there is no guarantee that the erroris distributed as a Gaussian, so the error term should be taken only as a roughindication of probable error.Suppose that you want to integrate a function g over a region W that is noteasy to sample randomly.
For example, W might have a very complicated shape.No problem. Just find a region V that includes W and that can easily be sampled(Figure 7.6.1), and then define f to be equal to g for points in W and equal to zerofor points outside of W (but still inside the sampled V ). You want to try to makeV enclose W as closely as possible, because the zero values of f will increase theerror estimate term of (7.6.1). And well they should: points chosen outside of Whave no information content, so the effective value of N , the number of points, isreduced. The error estimate in (7.6.1) takes this into account.General purpose routines for Monte Carlo integration are quite complicated(see §7.8), but a worked example will show the underlying simplicity of the method.Suppose that we want to find the weight and the position of the center of mass of anobject of complicated shape, namely the intersection of a torus with the edge of alarge box.
In particular let the object be defined by the three simultaneous conditions#($2z2 +x2 + y 2 − 3 ≤ 1(7.6.3)(torus centered on the origin with major radius = 4, minor radius = 2)x≥1y ≥ −3(7.6.4)(two faces of the box, see Figure 7.6.2). Suppose for the moment that the objecthas a constant density ρ.We want to estimate the following integrals over the interior of the complicatedobject:&&&&ρ dx dy dzxρ dx dy dzyρ dx dy dzzρ dx dy dz(7.6.5)The coordinates of the center of mass will be the ratio of the latter three integrals(linear moments) to the first one (the weight).In the following fragment, the region V , enclosing the piece-of-torus W , is therectangular box extending from 1 to 4 in x, −3 to 4 in y, and −1 to 1 in z.306Chapter 7.Random Numbersarea A∫fdxFigure 7.6.1.
Monte Carlo integration. Random points are chosen within the area A. The integral of thefunction f is estimated as the area of A multiplied by the fraction of random points that fall below thecurve f . Refinements on this procedure can improve the accuracy of the method; see text.y420124xFigure 7.6.2. Example of Monte Carlo integration (see text). The region of interest is a piece of a torus,bounded by the intersection of two planes.
The limits of integration of the region cannot easily be writtenin analytically closed form, so Monte Carlo is a useful technique.3077.6 Simple Monte Carlo Integration#include "nrutil.h"...n=...Set to the number of sample points desired.den=...Set to the constant value of the density.sw=swx=swy=swz=0.0;Zero the various sums to be accumulated.varw=varx=vary=varz=0.0;vol=3.0*7.0*2.0;Volume of the sampled region.for(j=1;j<=n;j++) {x=1.0+3.0*ran2(&idum);Pick a point randomly in the sampled rey=(-3.0)+7.0*ran2(&idum);gion.z=(-1.0)+2.0*ran2(&idum);if (z*z+SQR(sqrt(x*x+y*y)-3.0) < 1.0) {Is it in the torus?sw += den;If so, add to the various cumulants.swx += x*den;swy += y*den;swz += z*den;varw += SQR(den);varx += SQR(x*den);vary += SQR(y*den);varz += SQR(z*den);}}w=vol*sw/n;The values of the integrals (7.6.5),x=vol*swx/n;y=vol*swy/n;z=vol*swz/n;dw=vol*sqrt((varw/n-SQR(sw/n))/n);and their corresponding error estimates.dx=vol*sqrt((varx/n-SQR(swx/n))/n);dy=vol*sqrt((vary/n-SQR(swy/n))/n);dz=vol*sqrt((varz/n-SQR(swz/n))/n);A change of variable can often be extremely worthwhile in Monte Carlointegration.
Suppose, for example, that we want to evaluate the same integrals,but for a piece-of-torus whose density is a strong function of z, in fact varyingaccording toρ(x, y, z) = e5z(7.6.6)One way to do this is to put the statementden=exp(5.0*z);inside the if (...) block, just before den is first used.
This will work, but it isa poor way to proceed. Since (7.6.6) falls so rapidly to zero as z decreases (downto its lower limit −1), most sampled points contribute almost nothing to the sumof the weight or moments. These points are effectively wasted, almost as badlyas those that fall outside of the region W .
A change of variable, exactly as in thetransformation methods of §7.2, solves this problem. Letds = e5z dzso thats=1 5ze ,5z=1ln(5s)5(7.6.7)Then ρdz = ds, and the limits −1 < z < 1 become .00135 < s < 29.682. Theprogram fragment now looks like this308Chapter 7.Random Numbers#include "nrutil.h"...n=...Set to the number of sample points desired.sw=swx=swy=swz=0.0;varw=varx=vary=varz=0.0;ss=0.2*(exp(5.0)-exp(-5.0))Interval of s to be random sampled.vol=3.0*7.0*ssVolume in x,y,s-space.for(j=1;j<=n;j++) {x=1.0+3.0*ran2(&idum);y=(-3.0)+7.0*ran2(&idum);s=0.00135+ss*ran2(&idum);Pick a point in s.z=0.2*log(5.0*s);Equation (7.6.7).if (z*z+SQR(sqrt(x*x+y*y)-3.0) < 1.0) {sw += 1.0;Density is 1, since absorbed into definitionswx += x;of s.swy += y;swz += z;varw += 1.0;varx += x*x;vary += y*y;varz += z*z;}}w=vol*sw/n;The values of the integrals (7.6.5),x=vol*swx/n;y=vol*swy/n;z=vol*swz/n;dw=vol*sqrt((varw/n-SQR(sw/n))/n);and their corresponding error estimates.dx=vol*sqrt((varx/n-SQR(swx/n))/n);dy=vol*sqrt((vary/n-SQR(swy/n))/n);dz=vol*sqrt((varz/n-SQR(swz/n))/n);If you think for a minute, you will realize that equation (7.6.7) was useful onlybecause the part of the integrand that we wanted to eliminate (e5z ) was both integrableanalytically, and had an integral that could be analytically inverted.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
