#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "real.h"
#define rvector(n) (real*)calloc((n),sizeof(real))
#define ivector(n) (int*)calloc((n),sizeof(int))
#define tiny 0.000001

static int irec=0;

real strata(int n,real *a,real *b,real (*func)(real*),
	int m,real esq, void (*rnd)(int,real*,real*,real*))
{
int i,l; real fbar,f2bar,f,f2,big,*x,ssq,vol; x=rvector(n);

irec++; //if(irec>200){ printf("irec>200"); return 0;}

for(vol=1,i=0;i<n;i++) vol*=(b[i]-a[i]);
if (m<n) m=n; if (esq<tiny*tiny) esq=tiny*tiny;
for(fbar=f2bar=big=0, l=0, i=0;i<m;i++){
	(*rnd)(n,a,b,x); f=(*func)(x); fbar+=f; f2bar+=f*f;
	if(i<n){
		if(x[i]>(a[i]+b[i])/2) 	x[i]=x[i]-(b[i]-a[i])/2;
		else 				x[i]=x[i]+(b[i]-a[i])/2;
		f2=(*func)(x); fbar+=f2; f2bar+=f2*f2;
		if(fabs(f-f2)>big) {big=fabs(f-f2); l=i;}
		}
	}
fbar/=(m+n); f2bar/=(m+n); ssq=vol*vol*(f2bar-fbar*fbar)/(m+n-1);

if( ssq < esq )
	{
	 return vol*fbar;
	}
else{
	real al,bl,cl,temp;
	al=a[l]; bl=b[l]; cl=(al+bl)/2;
	m=m*0.707;
	if( ssq/esq >m ) m=(ssq/esq);
//	esq=esq*ssq/(ssq-esq);

	b[l]=cl;
	temp=strata(n,a,b,func,m,esq/2,rnd);
	b[l]=bl;

	a[l]=cl;
	temp+=strata(n,a,b,func,m,esq/2,rnd);
	a[l]=al;

	return (temp*ssq+vol*fbar*esq)/(ssq+esq);
	}
}