7.8 Adaptive and Recursive Monte Carlo Methods 321 void rebin(float rc,int nd,float r[],float xin[],float xi[]); static int i,it,j,k,mds,nd,ndo,ng,npg,ia[MXDIM+1],kg [MXDIM+1]; static float calls,dv2g,dxg,f,f2,f2b,fb,rc,ti,tsi,wgt,xjac,xn,xnd,xo; static float d[NDMX+1][MXDIM+1],di [NDMX+1][MXDIM+1],dt [MXDIM+1], dx [MXDIM+1],r [NDMX+1],x [MXDIM+1],xi [MXDIM+1][NDMX+1],xin [NDMX+1]; static double schi,si,swgt; Best make everything static,allowing restarts if (init <=0){ Normal entry.Enter here on a cold start. mds=ndo=1; Change to mds=0 to disable stratified sampling. for (j=1;j<=ndim;j++)xi[j][1]=1.0; i.e.,use importance sampling only. if (init <1)si=swgt=schi=0.0; Enter here to inherit the grid from a previous call,but not its answers. if (init <=2){ Enter here to inherit the previous grid and its nd=NDMX: answers. granted for 18881992 ng=1; if (mds){ Set up for stratification 1.800 (including this one) ng=(int)pow(ncall/2.0+0.25,1.0/ndim); mds=1; 1f(2*ng-lDM)>=0){ mds =-1; 7423 from NUMERICAL RECIPES IN npg=ng/NDMX+1; nd=ng/npg; ng=npg*nd; (North America to to users to make one paper for (k=1,i=1;i<-ndim;i++)k*=ng; by Cambridge University Press.Programs THE ART npg=IMAX(ncall/k,2); calls=(float)npg *(float)k; dxg=1.0/ng; for (dv2g=1,i=1;i<=ndim;i++)dv2g *dxg; dv2g=SQR(calls*dv2g)/npg/npg/(npg-1.0); xnd=nd; dxg *xnd; email xjac=1.0/calls; for (j=1;i<=ndim;j++){ dx[j]=regn[j+ndim]-regn[j]; xjac dx[j]; to directcustsen if(nd！=ndo)( Do binning if necessary for (i=1;i<=IMAX(nd,ndo);i++)r[i]=1.0; for (j=1;j<=ndim;j++)rebin(ndo/xnd,nd,r,xin,xi[j]); ndo=nd; @cambridge.org if (nprn >=0){ To order Numerical Recipes books or personaluse.Further reproduction,or 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN 0-521- printf("%s:ndim=%3d ncall=%8.0f\n", -431085 Input parameters for vegas",ndim,calls); printf("%28s it=%5d itmx=45d\n","",it,itmx); printf ("%28s nprn=%3d ALPH=%5.2f\n","",nprn,ALPH) printf("%28s mds=%3d nd=%4d\n","",mds,nd); Software. for (j=1;j<=ndim;j++){ printf("7%30s1[%2d]=%11.4gxu[%2d]=%11.4g\n", (outside North America) "j,regn[j],j,regn[j+ndim]) machine for (it=1:it<=itmx:it++){ Main iteration loop.Can enter here (init 3)to do an additional itmx iterations with all other parameters unchanged. t1=tsi=0.0; for (j=1;j<=ndim;j++) kg[j]=1; for(i=1;i<=nd;i+)d[i][j]=di[i][j]=0.0;7.8 Adaptive and Recursive Monte Carlo Methods 321 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). void rebin(float rc, int nd, float r[], float xin[], float xi[]); static int i,it,j,k,mds,nd,ndo,ng,npg,ia[MXDIM+1],kg[MXDIM+1]; static float calls,dv2g,dxg,f,f2,f2b,fb,rc,ti,tsi,wgt,xjac,xn,xnd,xo; static float d[NDMX+1][MXDIM+1],di[NDMX+1][MXDIM+1],dt[MXDIM+1], dx[MXDIM+1], r[NDMX+1],x[MXDIM+1],xi[MXDIM+1][NDMX+1],xin[NDMX+1]; static double schi,si,swgt; Best make everything static, allowing restarts. if (init <= 0) { Normal entry. Enter here on a cold start. mds=ndo=1; Change to mds=0 to disable stratified sampling, for (j=1;j<=ndim;j++) xi[j][1]=1.0; i.e., use importance sampling only. } if (init <= 1) si=swgt=schi=0.0; Enter here to inherit the grid from a previous call, but not its answers. if (init <= 2) { Enter here to inherit the previous grid and its nd=NDMX; answers. ng=1; if (mds) { Set up for stratification. ng=(int)pow(ncall/2.0+0.25,1.0/ndim); mds=1; if ((2*ng-NDMX) >= 0) { mds = -1; npg=ng/NDMX+1; nd=ng/npg; ng=npg*nd; } } for (k=1,i=1;i<=ndim;i++) k *= ng; npg=IMAX(ncall/k,2); calls=(float)npg * (float)k; dxg=1.0/ng; for (dv2g=1,i=1;i<=ndim;i++) dv2g *= dxg; dv2g=SQR(calls*dv2g)/npg/npg/(npg-1.0); xnd=nd; dxg *= xnd; xjac=1.0/calls; for (j=1;j<=ndim;j++) { dx[j]=regn[j+ndim]-regn[j]; xjac *= dx[j]; } if (nd != ndo) { Do binning if necessary. for (i=1;i<=IMAX(nd,ndo);i++) r[i]=1.0; for (j=1;j<=ndim;j++) rebin(ndo/xnd,nd,r,xin,xi[j]); ndo=nd; } if (nprn >= 0) { printf("%s: ndim= %3d ncall= %8.0f\n", " Input parameters for vegas",ndim,calls); printf("%28s it=%5d itmx=%5d\n"," ",it,itmx); printf("%28s nprn=%3d ALPH=%5.2f\n"," ",nprn,ALPH); printf("%28s mds=%3d nd=%4d\n"," ",mds,nd); for (j=1;j<=ndim;j++) { printf("%30s xl[%2d]= %11.4g xu[%2d]= %11.4g\n", " ",j,regn[j],j,regn[j+ndim]); } } } for (it=1;it<=itmx;it++) { Main iteration loop. Can enter here (init ≥ 3) to do an additional itmx iterations with all other parameters unchanged. ti=tsi=0.0; for (j=1;j<=ndim;j++) { kg[j]=1; for (i=1;i<=nd;i++) d[i][j]=di[i][j]=0.0;