package Lecture20;

class FuncC implements MathFunction {
    public double f(double x) {
        return x*x*x*x + 2;
    }
}

public class Simpson {
    public static double qsimp(MathFunction func, double a, double b) {
        double ost= -1.0E30;            // 
        double os= -1E30;
        for (int j=0; j < JMAX; j++) {
            double st= Trapezoid.trapzd(func, a, b, j+1);
            s= (4.0*st - ost)/3.0;      // Ver la ecuación Ref. Form. 4.2.4                 
            if (j > 4)                  // Evitar la primera convergencia espúrea
                if (Math.abs(s-os) < EPSILON*Math.abs(os) || (s==0.0 && os==0.0)) {
                    System.out.println("Iteraciones de Simpson: " + j);
                    return s; }
            os= s;
            ost= st;
        }
        System.out.println("Demasiados pasos en qsimp");
        return ERR_VAL;
    }
    private static double s;
    public static final double EPSILON= 1.0E-15;
    public static final int JMAX= 50;
    public static final double ERR_VAL= -1E10;

    public static void main(String[] args) {
        // Ejemplo sencillo solamente con trapezoide (ver RF, p. 137)                  
        System.out.println("Uso de trapezoide simple");
        int m= 20;   // Queremos 2^m+1 pasos
        int j= m+1;
        double ans= 0.0;
        for (j=0; j <=m; j++) {
            ans= Trapezoid.trapzd(new FuncC(), 0.0, 8.0, j+1);
            System.out.println("Iteración: " + (j+1) + " Integral: " + ans);
        }
        System.out.println("Integral: " + ans);
        
        System.out.println("Uso de Simpson");
        ans= qsimp(new FuncC(), 0.0, 8.0);
        System.out.println("Integral: " + ans);
        
        System.exit(0);        
    }
}

class Trapezoid {
    public static double trapzd(MathFunction func, double a, double b, int n) {
        if (n==1) {
            s= 0.5*(b-a)*(func.f(a)+func.f(b));
            return s; }
        else {
            int it= 1;
            for (int j= 0; j < n-2; j++)
                it *= 2;              // Subdivisiones
            double tnm= it;
            double delta= (b-a)/tnm;  // Espaciado de puntos que se van a agregar
            double x= a+0.5*delta;
            double sum= 0.0;
            for (int j= 0; j < it; j++) {
                sum += func.f(x);
                x+= delta; }
            s= 0.5*(s+(b-a)*sum/tnm);
            return s;
        }
    }
    private static double s;    // Conserva el cálculo actual de la integral. Sólo uno cada vez.
}