// CMPSC 201
// integration demo #1

#include <iostream>
#include <cmath>
using namespace std;

const double PI = 3.14159265;

double f(double x);
double trapezoid(double lowerLimit, double upperLimit, int n);
double simpson1_3(double lowerLimit, double upperLimit, int n);
double simpson3_8(double lowerLimit, double upperLimit, int n);

int main()
{
        double sum = 0.5;
	double sumT = 0.0;		// sum according to Trapezoid
	double sumS1_3 = 0.0;	// sum according to Simpson's 1/3
	double sumS3_8 = 0.0;	// sum according to Simpson's 3/8
	double h = 0.0;			// height
        int n = 0;				// number of segments
                
        cout << "  h  " << '\t' << " Trapezoid " << "\t" << " Simpsons(1/3) " << '\t' <<" Simpson's(3/8)" << endl;

        // use the property that the integral under entire f(x) is 1,
        // therefore the integral from -infinity to 0 is 0.5, so sum is initialized to 0.5 
        // and the integral is taken between 0 and 1.

        n = 1;
        h = 1.0;
        sumT = sum + trapezoid(0, 1, n);
        sumS1_3 = sum + simpson1_3(0, 1, n);
        sumS3_8 = sum + simpson3_8(0, 1, n);
        cout << h << '\t' << sumT<< '\t' << sumS1_3<< '\t' << sumS3_8 << endl;
        

        n = 10;
        h = 1.0/n;
        sumT = sum + trapezoid(0,1,n);
        sumS1_3 = sum + simpson1_3(0, 1, n);
        sumS3_8 = sum + simpson3_8(0, 1, n);
        cout << h << '\t' << sumT<< '\t' << sumS1_3<< '\t' << sumS3_8 << endl;
        

        n = 100;
        h = 1.0/n;
        sumT = sum + trapezoid(0,1,n);
        sumS1_3 = sum + simpson1_3(0, 1, n);
        sumS3_8 = sum + simpson3_8(0, 1, n);
        cout << h << '\t' << sumT<< '\t' << sumS1_3<< '\t' << sumS3_8 << endl;
        
        n = 1000;
        h = 1.0/n;
        sumT = sum + trapezoid(0, 1, n);
        sumS1_3 = sum + simpson1_3(0, 1, n);
        sumS3_8 = sum + simpson3_8(0, 1, n);
        cout << h << '\t' << sumT<< '\t' << sumS1_3<< '\t' << sumS3_8 << endl;

		return 0;
}// end of main

double trapezoid(double lowerLimit, double upperLimit, int n)
{
        double sum = 0.0;
	double area = 0.0;
	double a = 0.0;
	double b = 0.0;
	double h = 0.0;
        int i = 0; 

        h = (upperLimit - lowerLimit) / n;
        a = lowerLimit;
        b = lowerLimit + h;

        for (i = 1; i <= n; i++)
	{
		area = (b - a)/2 * (f(a) + f(b));
		sum = sum + area;
		a = b;
		b = a + h;
        }

        return sum;
}// end of trapezoid

double simpson1_3(double lowerLimit, double upperLimit, int n)
{
        double sum = 0.0;
	double area = 0.0;
	double a = 0.0;
	double b = 0.0;
	double h = 0.0;
        int i = 0; 

        h = (upperLimit - lowerLimit) / n;
        a = lowerLimit;
        b = lowerLimit + 2 * h;

        for (i = 1; i <= n; i += 2)
	{
		area = (b - a)/6 * (f(a) + 4 * f((b + a) / 2) + f(b));
		sum = sum + area;
		a = b;
		b = a + 2 * h;
        }

        return sum;
}// end of simpson1_3

double simpson3_8(double lowerLimit, double upperLimit, int n)
{
        double sum = 0.0;
	double area = 0.0;
	double a = 0.0;
	double b = 0.0;
	double h = 0.0;
        int i = 0; 

        h = (upperLimit - lowerLimit) / n;
        a = lowerLimit;
        b = lowerLimit + 3 * h;

        for (i = 1; i <= n; i += 3)
	{
		area = (b - a)/8 * (f(a) + 3 * f(a + h) + 3 * f(a+2 * h) + f(b));
		sum = sum + area;
		a = b;
		b = a + 3 * h;
        }

        return sum;
}// end of simpson3_8

double f(double x)
{
        double y = 0.0;
        
        y = 1 / sqrt(2 * PI) * exp(-(x * x) / 2);

        return y;
}// end of f