326 Chapter 7. Random Numbers unsigned long n,npre,nptl,nptr; int j,jb; float avel,varl; float fracl,fval; float rgl,rgm,rgr,s,sigl,siglb,sigr,sigrb; float sum,sumb,summ,summ2; float *fmaxl,*fmaxr,*fminl,*fminr; float *pt,*rmid; pt=vector(1,ndim); if (npts MNBS){ Too few points to bisect;do straight summ=summ2=0.0; Monte Carlo. http://www.nr. for (n=1;n<=npts;n++){ ranpt(pt,regn,ndim); fval=(*func)(pt); 83g summ fval: 188891992 summ2+=fval米fval; 1600 (including this one) *ave=summ/npts *var=FMAX(TINY,(summ2-summ*summ/npts)/(npts*npts)); else Do the preliminary (uniform)sampling. to any from NUMERICAL RECIPES IN rmid=vector(1,ndim); npre=LMAX((unsigned long)(npts*PFAC),MNPT); (Nort to make fmaxl=vector(1,ndim); fmaxr-vector(1,ndim); America server computer, one paper /Cambridge University Press. THE fminl=vector(1,ndim); fminr-vector(1,ndim); ART for (j=1;j<=ndim;j++) Initialize the left and right bounds for 1ran=(1ran*2661+36979)%175000; each dimension. s=SIGN(dith,(float)(iran-87500)); Programs rmid[j]=(0.5+s)*regn[j]+(0.5-s)*regn [ndim+j]; fminl [i]=fminr[i]=BIG; fmaxl[j]=fmaxr[j]=-BIG; for (n=1in<enpre;n++){ Loop over the points in the sample to dir ranpt(pt,regn,ndim); fval=(*func)(pt); for (j=1;j<=ndim;j++){ Find the left and right bounds for each ectcustser OF SCIENTIFIC COMPUTING(ISBN if (pt[j]<=rmid[j]) dimension. fminl[i]-FMIN(fminl[i],fval); fmaxl[j]=FMAX(fmaxl[j],fval); a 1988-1992 by Numerical Recipes 10-621 else fminr[i]=FMIN(fminr[i],fval); fmaxr[j]=FMAX(fmaxr [j],fval); -43108 (outside sumb=BIG; Choose which dimension jb to bisect jb=0: siglb=sigrb=1.0: North Software. for (j=1;j<=ndim;j++) if (fmaxl[j]fminl[j]&&fmaxr[j]fminr[j]){ sigl=FMAX(TINY,pow(fmaxl[j]-fminl[j],2.0/3.0)); sigr=FMAX(TINY,pow(fmaxr [j]-fminr [j],2.0/3.0)); sum=sigl+sigr; Equation (7.8.24),see text if (sum<=sumb){ sumb=sum; ib=i; siglb=sigl; sigrb=sigr;326 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). unsigned long n,npre,nptl,nptr; int j,jb; float avel,varl; float fracl,fval; float rgl,rgm,rgr,s,sigl,siglb,sigr,sigrb; float sum,sumb,summ,summ2; float *fmaxl,*fmaxr,*fminl,*fminr; float *pt,*rmid; pt=vector(1,ndim); if (npts < MNBS) { Too few points to bisect; do straight summ=summ2=0.0; Monte Carlo. for (n=1;n<=npts;n++) { ranpt(pt,regn,ndim); fval=(*func)(pt); summ += fval; summ2 += fval * fval; } *ave=summ/npts; *var=FMAX(TINY,(summ2-summ*summ/npts)/(npts*npts)); } else { Do the preliminary (uniform) sampling. rmid=vector(1,ndim); npre=LMAX((unsigned long)(npts*PFAC),MNPT); fmaxl=vector(1,ndim); fmaxr=vector(1,ndim); fminl=vector(1,ndim); fminr=vector(1,ndim); for (j=1;j<=ndim;j++) { Initialize the left and right bounds for iran=(iran*2661+36979) % 175000; each dimension. s=SIGN(dith,(float)(iran-87500)); rmid[j]=(0.5+s)*regn[j]+(0.5-s)*regn[ndim+j]; fminl[j]=fminr[j]=BIG; fmaxl[j]=fmaxr[j] = -BIG; } for (n=1;n<=npre;n++) { Loop over the points in the sample. ranpt(pt,regn,ndim); fval=(*func)(pt); for (j=1;j<=ndim;j++) { Find the left and right bounds for each if (pt[j]<=rmid[j]) { dimension. fminl[j]=FMIN(fminl[j],fval); fmaxl[j]=FMAX(fmaxl[j],fval); } else { fminr[j]=FMIN(fminr[j],fval); fmaxr[j]=FMAX(fmaxr[j],fval); } } } sumb=BIG; Choose which dimension jb to bisect. jb=0; siglb=sigrb=1.0; for (j=1;j<=ndim;j++) { if (fmaxl[j] > fminl[j] && fmaxr[j] > fminr[j]) { sigl=FMAX(TINY,pow(fmaxl[j]-fminl[j],2.0/3.0)); sigr=FMAX(TINY,pow(fmaxr[j]-fminr[j],2.0/3.0)); sum=sigl+sigr; Equation (7.8.24), see text. if (sum<=sumb) { sumb=sum; jb=j; siglb=sigl; sigrb=sigr; } }