subroutine plainmc(n,a,b,x,fun,m,res,acc,eps,point)
dimension a(n),b(n),x(n); external fun,point; parameter(nmax=1000000)

vol=1; do i=1,n; vol=vol*(b(i)-a(i)); end do

sum=0; sum2=0; do k=1,m
	call point(n,a,b,x); f=fun(n,x); sum=sum+f; sum2=sum2+f**2
end do

avg=sum/m; var=sum2/m-avg**2; err=vol*sqrt(var/m); res=vol*avg;
tol=acc+eps*res;

do while(err>tol)
	m=m+1; if(m>nmax) stop 'plainmc: m>nmax'
	call point(n,a,b,x); f=fun(n,x); sum=sum+f; sum2=sum2+f**2
	avg=sum/m; var=sum2/m-avg**2; err=vol*sqrt(var/m); res=vol*avg;
	tol=acc+eps*res;
end do

return; end