// Arup Guha
// 2/22/2012
// Trapezoid Rule

#include <stdio.h>
#include <math.h>

double trapezoid(double a, double b, int n);
double simpson(double a, double b, int n);
double f(double x);

int main() {
    double pi = 3.141592654;

    printf("%lf\n", trapezoid(-1,3,10000));
    printf("%lf\n", simpson(-1,3,10000));
    return 0;
}

// Adds up n trapezoids to approximate the area
// under the function f from [a,b].
double trapezoid(double a, double b, int n) {

    int i;

    // This is called h in the book.
    double step = (b - a)/n;

    double area = 0;

    // Add up each control point the
    // appropriate number of times.
    for (i=0; i<=n; i++) {

        double oldarea = area;
        // Ends are included once.
        if (i==0 || i == n)
            area += f(a+i*step);

        // Everything else is twice.
        else
            area += (2*f(a+i*step));

        printf("Just added %lf\n", area-oldarea);
    }

    // Now get the area of all the trapezoids...
    area = area*step/2;
    return area;

}

// Uses Simpson's Rule to approximate the area
// from [a,b] using n+1 sample points.
// Note: n must be even!!!
double simpson(double a, double b, int n) {

    int i;

    // This is called h in the book.
    double step = (b - a)/n;

    double area = 0;

    // Add up each control point the
    // appropriate number of times.
    for (i=0; i<=n; i++) {

        // Ends are included once.
        if (i==0 || i == n)
            area += f(a+i*step);

        // These are included five times.
        else if (i%2 == 1)
            area += (4*f(a+i*step));

        // Everything else is six times.
        else
            area += (2*f(a+i*step));
    }

    // Divide by 3 since we sample at three points.
    // Really, we would divide by 6, but we are adding two
    // regions of width step...
    area = area*step/3;
    return area;

}

// e^x for testing.
/*
double f(double x) {
    return exp(sin(x));
}
*/

double f(double x) {
    return exp(x*x);
}