using System;
using static System.Console;
using static System.Math;
using static System.Double;

public partial class montecarlo{

public static double[] strata
(Func<vector,double> f, double[] a, double[] b,
double acc=1e-2, double eps=1e-2, Random RND=null,
double old_N=0,double old_average=0)
{
if(RND==null)RND=new Random(1);
int dim=a.Length;
int N=16*dim;

double volume=1; for(int i=0;i<dim;i++) volume*=(b[i]-a[i]);

var x =new vector(dim);

if(old_N==0){
	old_N=N;
	old_average=0;
	for(int n=0;n<N;n++){
		for(int i=0;i<x.size;i++) x[i]=a[i]+RND.NextDouble()*(b[i]-a[i]);
		double fx=f(x);
		old_average+=fx;
		}
	old_average/=N;
	}

var aL=new vector(dim);
var aR=new vector(dim);
var vL=new vector(dim);
var vR=new vector(dim);
var nL=new int[dim];
var nR=new int[dim];
var sum_f_L=new vector(dim);
var sum_f2_L=new vector(dim);
var sum_f_R=new vector(dim);
var sum_f2_R=new vector(dim);

for(int n=0;n<N;n++){
	for(int i=0;i<x.size;i++) x[i]=a[i]+RND.NextDouble()*(b[i]-a[i]);
	double fx=f(x);
	for(int i=0;i<dim;i++){ // sub-sums
		if(x[i]<(a[i]+b[i])/2){nL[i]++;sum_f_L[i]+=fx;sum_f2_L[i]+=fx*fx;}
		else                  {nR[i]++;sum_f_R[i]+=fx;sum_f2_R[i]+=fx*fx;}
		}
	}

double vmax=0;int i_bisect=RND.Next(0,dim);
for(int i=0;i<dim;i++){
	if(nL[i]>0){
		aL[i]=sum_f_L[i]/nL[i];
		vL[i]=sum_f2_L[i]/nL[i]-Pow(aL[i],2);
		}
	if(nR[i]>0){
		aR[i]=sum_f_R[i]/nR[i];
		vR[i]=sum_f2_R[i]/nR[i]-Pow(aR[i],2);
		}
	double measure=Abs(aL[i]-aR[i]);
	//double measure=vL[i]+vR[i];
	//double measure=Max(vL[i],vR[i]);
	if(measure>vmax){vmax=measure;i_bisect=i;}
	}

double average=(sum_f_L[i_bisect]+sum_f_R[i_bisect])/N;
double integ=(average*N+old_average*old_N)/(N+old_N)*volume;
double error=Abs(average-old_average)*volume;

if(error<acc+Abs(integ)*eps)return new double[] {integ,error};

var aa = (double[])a.Clone();aa[i_bisect]=(a[i_bisect]+b[i_bisect])/2;
var bb = (double[])b.Clone();bb[i_bisect]=(a[i_bisect]+b[i_bisect])/2;

double[] result_L =
	strata(f,a,bb,acc/Sqrt(2.0),eps,RND,nL[i_bisect],aL[i_bisect]);
double[] result_R =
	strata(f,aa,b,acc/Sqrt(2.0),eps,RND,nR[i_bisect],aR[i_bisect]);

integ=result_L[0]+result_R[0];
error=Sqrt(Pow(result_L[1],2)+Pow(result_R[1],2));

return new double[] {integ,error};

}//miser
}//class