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

public static class strata{

public static int Nmin=42;
public static Random rnd;

public static double[] integrate(
Func<vector,double> f,vector a, vector b, int N)
	{
	rnd=new Random();
	Nmin=10*a.size;
	double V=1;
	for(int i=0;i<a.size;i++) V*=(b[i]-a[i]);
	double[] r=sample(f,a,b,N);
	double integ=r[0]*V;
	double error=r[1]*V/Sqrt(N);
	return new double[] {integ,error};
	}
	

public static double[] sample(
	Func<vector,double> f,vector a, vector b, int N){

Error.Write($"sample: N={N} Nmin={Nmin}\n");
if(N<=0)return new double[] {0,0};

int dim=a.size;
var x =new vector(dim);
Action<vector> randomx=(vector v)=>{
	for(int i=0;i<dim;i++)v[i]=a[i]+rnd.NextDouble()*(b[i]-a[i]);
	};

if(N<=Nmin) {
        double sum=0,sum2=0;
        for(int i=0;i<N;i++){
                randomx(x);
                double fx=f(x);
                sum+=fx;
                sum2+=fx*fx;
                }
        double average=sum/N, variance=sum2/N-average*average;
        return new double[] {average,variance};
	}

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];

int n=Max(Nmin,N/10);
for(int k=0;k<n;k++){
	randomx(x);
	double fx=f(x);
	for(int i=0;i<dim;i++){
		if(x[i]<(a[i]+b[i])/2){nL[i]++;aL[i]+=fx;vL[i]+=fx*fx;}
		else                  {nR[i]++;aR[i]+=fx;vR[i]+=fx*fx;}
		}
	}

double vmax=0;int imax=0;
for(int i=0;i<dim;i++){
	if(nL[i]>0){
		aL[i]/=nL[i];
		vL[i]=vL[i]/nL[i]-Pow(aL[i],2);
		}
	if(nR[i]>0){
		aR[i]/=nR[i];
		vR[i]=vR[i]/nR[i]-Pow(aR[i],2);
		}
	double measure=Abs(aL[i]-aR[i]);
	if(measure>vmax){vmax=measure;imax=i;}
	}

var aa=a.copy();aa[imax]=(a[imax]+b[imax])/2;
var bb=b.copy();bb[imax]=(a[imax]+b[imax])/2;

int NL,NR;
if(vL[imax]+vR[imax]==0) NL=(N-n)/2;
else NL=(int)Round((N-n)*vL[imax]/(vL[imax]+vR[imax]));
NR=N-n-NL;

double[] rL=sample(f,a,bb,NL);
double[] rR=sample(f,aa,b,NR);
Error.Write($"rL = {(float)rL[0]} {(float)rL[1]}\n");
Error.Write($"oL = {(float)aL[imax]} {(float)vL[imax]}\n");
Error.Write($"rR = {(float)rR[0]} {(float)rR[1]}\n");
Error.Write($"oR = {(float)aR[imax]} {(float)vR[imax]}\n");
Error.Write($"NL+nL= {NL+nL[imax]}\n");
Error.Write($"NR+nR= {NR+nR[imax]}\n");
double AL=(NL*rL[0]+nL[imax]*aL[imax])/(NL+nL[imax]);
double AR=(NR*rR[0]+nR[imax]*aR[imax])/(NR+nR[imax]);
double VL=(NL*rL[1]+nL[imax]*vL[imax])/(NL+nL[imax]);
double VR=(NR*rR[1]+nR[imax]*vR[imax])/(NR+nR[imax]);
Error.Write($" L = {(float)AL} {(float)VL} ");
Error.Write($" L = {(float)AR} {(float)VR}\n");
return new double[] {0.5*AL+0.5*AR,0.5*VL+0.5*VR};
}//integrate
}//strata