/* This class has two methods for approximating the integral of a function
   over a specified range, one iterative (iterIntegrate()) and one recursive
   (recIntegrate()).  Each of them computes the integral of any given
   one-dimensional function.  The main() method tests these methods on a few
   examples.

   To represent a given function, like x^2, we define a class whose yVal(x)
   method returns x^2 (class Square below).  All such classes must implement
   the Function interface defined below.  See the code and comments for these
   classes to see what this means and how it's done.

   Written by Fatima Ahmad in May 1999; tinkered with by Jacob Eliosoff in
   October 1999.  */
class Integrate
{

    /* Returns an (iteratively computed) approximation of the integral of f
       over the range (xMin, xMax).  The smaller the step size, the more
       accurate the estimate.

       How it works:  The integral of f from xMin to xMax is the area under f
       (that is, between f and the x-axis) between those values.  We can
       approximate this area by dividing it into many narrow vertical strips.
       The step size determines the width of each strip, and if it's small
       enough, each strip looks like a rectangle.  The area of a rectangle is
       easy to compute, so we can approximate the total area by adding up the
       areas of all the rectangles.  */
    public static double iterIntegrate(Function f,
                                       double xMin, double xMax, double step)
    {
        double result, x, xMiddle, stripArea;

        result = 0.0;
        for (x = xMin; x < xMax; x += step) {
            /* The height of f at the current x-value is f.yVal(x); at x+step,
               the height is f.yVal(x+step).  For small values of step, these
               two y-values will be very close.  So we can approximate the area
               of the strip from x to x+step by calculating the area of the
               rectangle with about the same height, say the value of f at the
               midpoint between x and x+step:  */
            xMiddle = x + (step/2);
            stripArea = step * f.yVal(xMiddle);

            /* Having done that, we add the rectangle's area to the total so
               far:  */
            result += stripArea;
        }

        return result;
    }

    /* Returns a (recursively computed) approximation of the integral of f over
       the range (xMin, xMax).  The smaller the value of epsilon, the more
       accurate the estimate.

       How it works:  Again, we want to compute the area under f between xMin
       and xMax.  This time, rather than dividing into many narrow strips, we
       compute the value of f at xMin, xMax, and the midpoint xMiddle.  If
       these three points form a nearly straight line, then we conclude that
       over this range f can be approximated by a straight line, allowing us to
       approximate its area as that of the trapezoid under the straight line
       from (xMin, yMin) to (xMax, yMax).  (This is the base case.)  If xMiddle
       doesn't lie within epsilon of this line, we conclude that f is too
       erratic over this range to be approximated by a trapezoid, so we divide
       further, splitting the range into a left half and a right half and
       approximating their areas recursively.

       The point of this algorithm is that its approximation uses finer strips
       where the function is erratic, and wider ones where the function is
       smooth, whereas the iterative algorithm above uses strips of the same
       width over the entire range.  */
    public static double recIntegrate(Function f,
                                      double xMin, double xMax, double epsilon)
    {
        double minVal, maxVal, midVal, midEst, xMiddle;
        
        xMiddle = (xMin + xMax)/2;
        
        minVal = f.yVal(xMin);
        maxVal = f.yVal(xMax);
        midVal = f.yVal(xMiddle);
        
        midEst = minVal + (maxVal-minVal)/2; //find point on estimation line
        
        //check against epsilon:
        if (Math.abs(midVal-midEst) <= epsilon) {
            return .5*(minVal + maxVal)*(xMax - xMin);
        }
        else {
            return (recIntegrate(f, xMin, xMiddle, epsilon) + recIntegrate(f, xMiddle, xMax, epsilon));
        }
    }

    /* Tests the two integration methods above on the Function objects defined
       below.  */
    public static void main(String args[])
    {
        /* Setting up 8 test calls to the integration methods above.  Because
           Identity, Square, etc implement the Function interface, objects of
           these types are also considered Functions, and can be stored in
           arrays expecting elements of type Function, assigned to Function
           variables, and so on.  See comments below for the Function
           interface.  */
        Function[] functions = {
            new Identity(), new Square(),   new Cube(),     new F1(),
            new F2(),       new F2(),       new Sine(),     new Gaussian(),
        };
        double[] xMins = {
            0.0,    -1.0,   0.0,    -2.5,   -2.5,   -2.5,   0,    -3.0,
        };
        double[] xMaxes = {
            1.0,    1.0,    1.0,    2.5,    2.5,    2.5,    Math.PI, 3.0,
        };
        double[] steps = {
            0.01,   0.01,   0.01,   0.01,   0.01,   0.1,    0.01,   0.01,
        };
        double[] epsilons = {
            0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.01,   0.0001, 0.0001,
        };
        Function f;
        double xMin, xMax, step, epsilon;

        for (int i = 0; i < functions.length; ++i) {
            f = functions[i];
            xMin = xMins[i];
            xMax = xMaxes[i];

            /* String-concatenating f like this will cause f.toString() to be
               called (see toString() methods below):  */
            System.out.println("Integral of  " + f +
                               "  from " + xMin + " to " + xMax + ":");
            step = steps[i];
            System.out.println("    Iterative method (step = " + step +
                               "):  " + iterIntegrate(f, xMin, xMax, step));
            epsilon = epsilons[i];
            System.out.println("    Recursive method (epsilon = " + epsilon +
                               "):  " + recIntegrate(f, xMin, xMax, epsilon));
        }
    }

}


/* This is the interface which a class implements if it wants to be treated as
   a Function.  For example, Identity (below) implements Function, meaning that
   it provides bodies for all methods declared in the Function interface.
   (This particular interface declares only one method, but if more were
   declared, Identity and other Function implementers would have to provide
   bodies for all of them.)

   Interfaces like this allow us to define a general type which can refer to
   any of several different classes.  So, we can write a method which takes an
   object of type Function and call its yVal() method, without needing to
   know exactly which class it belongs to.  The integration methods above are
   like this.  We can define new classes implementing Function and pass them to
   iterIntegrate().  Without interfaces, it would be hard to write the
   integrate methods without hardcoding which functions it could integrate.  */
interface Function
{

    /* The one method declared by this interface.  Every class which claims at
       the top that it "implements Function" must provide a definition (body)
       for this method, returning the appropriate y-value for the given x.  */
    public double yVal(double x);

}


/* A very simple Function class representing the identity function.  */
class Identity
    implements Function
{

    /* Returns x.  */
    public double yVal(double x)
    {
        return x;
    }

    /* Overriding the special toString() method, in any class, determines what
       String will be substituted when an object of that class is converted to
       a String (eg, by being printed out by System.out.println()).  We
       override it here so that the String printed out in Integrate.main()
       above will properly describe the particular Function f in use.  */
    public String toString()
    {
        return "x";
    }

}


/* A Function class representing the x^2 function.  */
class Square
    implements Function
{

    /* Returns the appropriate y-value for the given x-value.  See above.  */
    public double yVal(double x)
    {
        return (x * x);
    }

    /* Returns a String describing this function.  See above.  */
    public String toString()
    {
        return "x^2";
    }

}


/* A Function class representing the x^3 function.  */
class Cube
    implements Function
{
 
    /* Returns the appropriate y-value for the given x-value.  See above.  */
   public double yVal(double x)
    {
        return (x * x * x);
    }

    /* Returns a String describing this function.  See above.  */
    public String toString()
    {
        return "x^3";
    }

}


/* A Function class representing the quartic polynomial shown below.  */
class F1
    implements Function
{

    /* Returns the appropriate y-value for the given x-value.  See above.  */
    public double yVal(double x)
    {
        return (Math.pow(x, 4) - (2.0 * Math.pow(x, 2)));
    }

    /* Returns a String describing this function.  See above.  */
    public String toString()
    {
        return "x^4 - 2x^2";
    }

}


/* A Function class representing the quintic polynomial shown below.  */
class F2
    implements Function
{

    /* Returns the appropriate y-value for the given x-value.  See above.  */
    public double yVal(double x)
    {
        return (Math.pow(x, 5) - Math.pow(x, 4) - (6.0 * Math.pow(x, 3)) +
                (4.0 * Math.pow(x, 2)) + (8.0 * x));
    }

    /* Returns a String describing this function.  See above.  */
    public String toString()
    {
        return "x^5 - x^4 - 6x^3 + 4x^2 + 8x";
    }

}


/* A Function class representing the sine function.  */
class Sine
    implements Function
{

    /* Returns the appropriate y-value for the given x-value.  See above.  */
    public double yVal(double x)
    {
        return Math.sin(x);
    }

    /* Returns a String describing this function.  See above.  */
    public String toString()
    {
        return "sin(x)";
    }

}


/* A Function class representing the Gaussian function, e^(-x^2).  */
class Gaussian
    implements Function
{

    /* Returns the appropriate y-value for the given x-value.  See above.  */
    public double yVal(double x)
    {
        return Math.exp(-(x * x));
    }

    /* Returns a String describing this function.  See above.  */
    public String toString()
    {
        return "e^(-x^2)";
    }

}