// CMPSC 201
// integration demo #2


#include <iostream>
using namespace std;

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 sumT = 0.0;		// integral according to Trapezoid
	double sumS1_3 = 0.0;		// integral according to Simpson's 1/3
	double sumS3_8 = 0.0;		// integral according to Simpson's 3/8
	double a = 0.0;			// lower limit
	double b = 0.0;			// upper limit
        int n = 0;			// number of segments

	cout << "The function to integrate is f(x) = 1/x" << endl;
	cout << "Enter lower limit (a): ";
	cin >> a;
	cout << "Enter upper limit (b): ";
	cin >> b;
	cout << "Enter odd number of segments (n, will compute with n+1 for Simpson's 1/3): ";
	cin >> n;
                
        cout << '\t' << " Trapezoid " << "\t" << " Simpsons(1/3) " << '\t' <<" Simpson's(3/8)" << endl;

        sumT = trapezoid(a, b, n);
        sumS1_3 = simpson1_3(a, b, n);
        sumS3_8 = simpson3_8(a, b, n);
        cout << '\t' << sumT << '\t' << sumS1_3 << '\t' << sumS3_8 << endl;

	return 0;
}// end of main

double trapezoid(double lowerLimit, double upperLimit, int n)
{						// a different trapezoid function from demo #1
        double sum = 0.0;
		double h = 0.0;
        int i = 0; 

        h = (upperLimit - lowerLimit) / n;
        sum = f(lowerLimit) + f(upperLimit);

        for (i = 1; i < n; i++)
	{
		sum += 2 * f(lowerLimit + i * h);
        }

        return h / 2.0 * sum;
}// end of trapezoid

double simpson1_3(double lowerLimit, double upperLimit, int n)
{
   // precondition: n is even

	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)
{
	// precondition: n is odd

	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 / x;

        return y;
}// end of f