7.7 Quasi-(that is,Sub-)Random Sequences 313 in=0; if (iv[1]!1)return; fac=1.0/(1L <MAXBIT) for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM)iu[j]&iv[k]; To allow both 1D and 2D addressing. for (k=1;k<=MAXDIM;k++){ for (j=1;j<=mdeg[k];j++)iu[j][k]<<(MAXBIT-j); Stored values only require normalization. for (j=mdeg[k]+1;j<=MAXBIT;j++) Use the recurrence to get other val- ipp=ip[k]; ues. i=iu[j-mdeg[k]][k]; i =(i>>mdeg[k]); for(1=mdeg[k]-1;1>=1;1--){ if (ipp 1)i iu[j-1][k] 1pp>=1; 鱼 18881892 iu[j][k]=i; 2 g NUMERICAL else Calculate the next vector in the se- imsin++; quence. for (j=1;j<=MAXBIT;j++){ Find the rightmost zero bit. RECIPES if (!(im&1))break; 1m>>=1; 令 if (j MAXBIT)nrerror("MAXBIT too small in sobseg"); im=(j-1)*MAXDIM: Americ computer Press. for (k=1;k<=IMIN(*n,MAXDIM);k++){ XOR the appropriate direction num- ix[k]“=1v[im+k]; ber into each component of the x[k]=ix[k]*fac; vector and convert to a floating number. SCIENTIFIC How good is a Sobol'sequence,anyway?For Monte Carlo integration of a smooth 6 function in n dimensions,the answer is that the fractional error will decrease with N,the number of samples,as(In N)"/N,i.e.,almost as fast as 1/N.As an example,let us integrate a function that is nonzero inside a torus (doughnut)in three-dimensional space.If the major radius of the torus is Ro,the minor radial coordinate r is defined by 1920 COMPUTING (ISBN r=(【e2+-P+2e (7.7.4) 10621 Let us try the function Numerica uctio 431086 f(x,y,z) cos r<To Recipes (7.7.5) 0 r20 (outside which can be integrated analytically in cylindrical coordinates,giving North Software. dr dy dz f(,y,z)=2n2a2Ro (7.7.6) With parameters Ro =0.6,ro =0.3,we did 100 successive Monte Carlo integrations of equation (7.7.4),sampling uniformly in the region-1 <y,z<1,for the two cases of uncorrelated random points and the Sobol'sequence generated by the routine sobseq.Figure 7.7.2 shows the results,plotting the r.m.s.average error of the 100 integrations as a function of the number of points sampled.(For any single integration,the error of course wanders from positive to negative,or vice versa,so a logarithmic plot of fractional error is not very informative.)The thin,dashed curve corresponds to uncorrelated random points and shows the familiar N-1/2 asymptotics.The thin,solid gray curve shows the result for the Sobol'7.7 Quasi- (that is, Sub-) Random Sequences 313 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). in=0; if (iv[1] != 1) return; fac=1.0/(1L << MAXBIT); for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM) iu[j] = &iv[k]; To allowboth 1D and 2D addressing. for (k=1;k<=MAXDIM;k++) { for (j=1;j<=mdeg[k];j++) iu[j][k] <<= (MAXBIT-j); Stored values only require normalization. for (j=mdeg[k]+1;j<=MAXBIT;j++) { Use the recurrence to get other val￾ipp=ip[k]; ues. i=iu[j-mdeg[k]][k]; i ^= (i >> mdeg[k]); for (l=mdeg[k]-1;l>=1;l--) { if (ipp & 1) i ^= iu[j-l][k]; ipp >>= 1; } iu[j][k]=i; } } } else { Calculate the next vector in the se￾im=in++; quence. for (j=1;j<=MAXBIT;j++) { Find the rightmost zero bit. if (!(im & 1)) break; im >>= 1; } if (j > MAXBIT) nrerror("MAXBIT too small in sobseq"); im=(j-1)*MAXDIM; for (k=1;k<=IMIN(*n,MAXDIM);k++) { XOR the appropriate direction num￾ber into each component of the vector and convert to a floating number. ix[k] ^= iv[im+k]; x[k]=ix[k]*fac; } } } How good is a Sobol’ sequence, anyway? For Monte Carlo integration of a smooth function in n dimensions, the answer is that the fractional error will decrease with N, the number of samples, as (ln N) n/N, i.e., almost as fast as 1/N. As an example, let us integrate a function that is nonzero inside a torus (doughnut) in three-dimensional space. If the major radius of the torus is R0, the minor radial coordinate r is defined by r =
[(x2 + y2 ) 1/2 − R0] 2 + z2 1/2 (7.7.4) Let us try the function f(x, y, z) =    1 + cos πr2 a2  r<r0 0 r ≥ r0 (7.7.5) which can be integrated analytically in cylindrical coordinates, giving  dx dy dz f(x, y, z)=2π2 a2 R0 (7.7.6) With parameters R0 = 0.6, r0 = 0.3, we did 100 successive Monte Carlo integrations of equation (7.7.4), sampling uniformly in the region −1 < x, y, z < 1, for the two cases of uncorrelated random points and the Sobol’ sequence generated by the routine sobseq. Figure 7.7.2 shows the results, plotting the r.m.s. average error of the 100 integrations as a function of the number of points sampled. (For any single integration, the error of course wanders from positive to negative, or vice versa, so a logarithmic plot of fractional error is not very informative.) The thin, dashed curve corresponds to uncorrelated random points and shows the familiar N−1/2 asymptotics. The thin, solid gray curve shows the result for the Sobol’