#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

double f(double x);
double F(double x);
double trapInt();
double rombInt();
double intPow(double x, int n);

double a = 10;
double b = 15;

int main() {
	
	cout << scientific << showpoint << setprecision(12)
			<< "Welcome to the Show!\n\n";
	
	cout << "Symbolic:\t" << F(b) - F(a) << endl;
	
	cout << "Trapezoid:\t" << trapInt() << endl;
	cout << "Romberg:\t" << rombInt() << endl;

	return 0;
}

/* returns f(x) for given x */
double f(double x) {
	return 2.1*log(x) - 3.0*sin(2*x) - 2.7*cos(4.5*x) - 0.3;
}

double F(double x) {
	return 2.1*(x*log(x)-x) + 1.5*cos(2*x) - 0.6*sin(4.5*x) - 0.3*x;
}

double trapInt() {
	
	double h, sum, x;
	int n;
	
	/* symbolic calculation gives that f''(q)<=58.994808 for 10<q<15 */
	h = sqrt(0.5e-6 / (5./12. * 58.994808));
	n = (int) ceil( (b-a) / h );
	
//	cout << "trap intervals: " << n << endl;
	
	sum = 0.5 * (f(a)+f(b));
	
	for (int i = 1; i <= n-1; i++) {
		x = a + i*h;
		sum += f(x);
	}
	
	sum *= h;
	
	return sum;
}

double rombInt() {
	
	int maxRows = 20;
	double epsilon = 0.5e-6;
	double R[maxRows][maxRows];
	double sum;
	
	double h = b - a;
	R[0][0] = h/2*(f(a)+f(b));
//	cout << R[0][0] << endl;

	for (int n = 1; n <= maxRows; n++) {
		h /= 2;
		sum = 0;
		
		for (int i = 1; i <= intPow(2,n)-1; i+=2) {
			sum += f(a+i*h);
		}
		
		R[n][0] = 0.5*R[n-1][0]+h*sum;
		
//		cout << "n = " << n << "\tR[n][0] = " << R[n][0] << endl;
		for (int m = 1; m <= n; m++) {
//			if ((R[n][m-1]-R[n-1][m-1]) < epsilon)
//				return R[n][m-1];
			R[n][m] = R[n][m-1];
//				+ (R[n][m-1]-R[n-1][m-1])/(intPow(4,m)-1);
		}
	}
	
	return 999;
}

double intPow(double x, int n) {
	if (n > 1)
		return x*intPow(x, n-1);
	else
		return x;
}