//======================================================
// Numerical Analysis package in java
// Adaptive simpson integration
// error control :
// 1.0/15.0*Math.abs(I1-is)<eps 
// Dr. Turhan Coban
// =====================================================
import java.io.*;

class f1 extends f_x
{
  double func(double x)
  {return x*x;}
}


class integ10
{
public static double adaptive_simpson_integral(f_x ff,double a,double b)
{return adaptive_simpson_integral(ff,a,b,1.0e-20,30);}

public static double adaptive_simpson_integral(f_x ff,double a,double b,double eps)
{return adaptive_simpson_integral(ff,a,b,eps,30);}

public static double adaptive_simpson_integral(f_x ff,double a,double b,int MAX_ITERATION)
{return adaptive_simpson_integral(ff,a,b,1.0e-12,MAX_ITERATION);}

public static double adaptive_simpson_integral(f_x ff,double a,double b,double eps,int MAX_ITERATION)
{   int i;
    double h=(b-a);
    double h2=h/4.0;
	double x[]=new double[5];
    double f[]=new double[5];
    for(i=0;i<5;i+=2) {x[i]=a+i*h2;f[i]=ff.func(x[i]);}
	double sum=h*(f[0]+4.0*f[2]+f[4])/6.0; 
	double Q=0;
    Q=simpson1_3(ff,x,f,sum,eps,0,MAX_ITERATION,Q);
    return Q;	
}
	
public static double simpson1_3(f_x ff,double x[],double f[],double sum,double eps,int k,int MAX_ITERATION,double Q)
{       //Newton-cotes integral 3 points
        //adaptive Simpson 1/3 rule
        double h=(x[4]-x[0]);
        double h1=h/2.0;
        double h2=h/4.0;        
        int p,i;
        double x1[]=new double[5];
        double f1[]=new double[5];
        double x2[]=new double[5];
        double f2[]=new double[5];
        k++;
        for(i=1;i<5;i+=2) {x[i]=x[0]+i*h2;f[i]=ff.func(x[i]);}
        double sum1=h1*(f[0]+4.0*f[1]+f[2])/6.0;
        double sum2=h1*(f[2]+4.0*f[3]+f[4])/6.0;     
        x1[0]=x[0];f1[0]=f[0];
        x1[2]=x[1];f1[2]=f[1];
        x1[4]=x[2];f1[4]=f[2];
        x2[0]=x[2];f2[0]=f[2];
        x2[2]=x[3];f2[2]=f[3];
        x2[4]=x[4];f2[4]=f[4];      
        if(1.0/15.0*Math.abs(sum1+sum2-sum)<eps || k>MAX_ITERATION) Q+=sum1+sum2;
        else Q+=simpson1_3(ff,x1,f1,sum1,eps,k,MAX_ITERATION,Q)+simpson1_3(ff,x2,f2,sum2,eps,k,MAX_ITERATION,Q);
        return Q;
}

  public static void main(String args[]) throws IOException
  {
  f1 b1=new f1();
  System.out.println("integral of class         f1 : "+adaptive_simpson_integral(b1,0.0,1.0));
}
}
