#include <stdlib.h>#include <math.h>#define NRANSI#include "nrutil.h"#define PFAC 0.1#define MNPT 15#define MNBS 60#define TINY 1.0e-30#define BIG 1.0e30static long iran=0;void miser(float (*func)(float []), float regn[], int ndim, unsigned long npts,	float dith, float *ave, float *var){	void ranpt(float pt[], float regn[], int n);	float *regn_temp;	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) {		summ=summ2=0.0;		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 {		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++) {			iran=(iran*2661+36979) % 175000;			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++) {			ranpt(pt,regn,ndim);			fval=(*func)(pt);			for (j=1;j<=ndim;j++) {				if (pt[j]<=rmid[j]) {					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;		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;				if (sum<=sumb) {					sumb=sum;					jb=j;					siglb=sigl;					sigrb=sigr;				}			}		}		free_vector(fminr,1,ndim);		free_vector(fminl,1,ndim);		free_vector(fmaxr,1,ndim);		free_vector(fmaxl,1,ndim);		if (!jb) jb=1+(ndim*iran)/175000;		rgl=regn[jb];		rgm=rmid[jb];		rgr=regn[ndim+jb];		fracl=fabs((rgm-rgl)/(rgr-rgl));		nptl=(unsigned long)(MNPT+(npts-npre-2*MNPT)*fracl*siglb			/(fracl*siglb+(1.0-fracl)*sigrb));		nptr=npts-npre-nptl;		regn_temp=vector(1,2*ndim);		for (j=1;j<=ndim;j++) {			regn_temp[j]=regn[j];			regn_temp[ndim+j]=regn[ndim+j];		}		regn_temp[ndim+jb]=rmid[jb];		miser(func,regn_temp,ndim,nptl,dith,&avel,&varl);		regn_temp[jb]=rmid[jb];		regn_temp[ndim+jb]=regn[ndim+jb];		miser(func,regn_temp,ndim,nptr,dith,ave,var);		free_vector(regn_temp,1,2*ndim);		*ave=fracl*avel+(1-fracl)*(*ave);		*var=fracl*fracl*varl+(1-fracl)*(1-fracl)*(*var);		free_vector(rmid,1,ndim);	}	free_vector(pt,1,ndim);}#undef MNPT#undef MNBS#undef TINY#undef BIG#undef PFAC#undef NRANSI