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

using namespace std;

//求積：梯形公式
double Func_Integral_Trapezoid
(double lo, double hi, double(*Func)(double), int n = 1000)
{//lo:下限；hi:上限；Func:函式；n:等分數

	if (n <= 0) n = 100;

	double x;
	double step = (hi - lo) / n;

	double result = 0.0;
	x = lo;
	for (int i = 1; i < n; i++) {
		x += step;
		result += Func(x);
	}
	result += (Func(lo) + Func(hi)) / 2;
	result *= step;

	return result;
}

//求積：Simpson公式
double Func_Integral_Simpson
(double lo, double hi, double(*Func)(double), int n = 1000)
{//lo:下限；hi:上限；Func:函式；n:等分數

	if (n <= 0) n = 100;

	double x;
	double step = (hi - lo) / n;

	double result1 = 0.0;
	x = lo;
	for (int i = 1; i < n; i++) {
		x += step;
		result1 += Func(x);
	}
	result1 *= 2;

	double result2 = 0.0;
	x = lo + step / 2;
	for (int i = 0; i < n; i++) {
		result2 += Func(x);
		x += step;
	}
	result2 *= 4;

	double result = result1 + result2 + Func(lo) + Func(hi);
	result *= step / 6;

	return result;
}

//求積：Cotes公式
double Func_Integral_Cotes
(double lo, double hi, double(*Func)(double), int n = 1000)
{//lo:下限；hi:上限；Func:函式；n:等分數

	if (n <= 0) n = 100;

	double x;
	double step = (hi - lo) / n;

	double result1 = 0.0;
	x = lo;
	for (int i = 1; i < n; i++) {
		x += step;
		result1 += Func(x);
	}
	result1 *= 14;

	double result2 = 0.0;
	x = lo + step / 2;

	double result3 = 0.0;
	double x1 = lo + step / 4;
	double x2 = lo + step / 4 * 3;

	for (int i = 0; i < n; i++) {
		result2 += Func(x);
		result3 += Func(x1) + Func(x2);
		x += step;
		x1 += step;
		x2 += step;
	}
	result2 *= 12;
	result3 *= 32;

	double result4 = (Func(lo) + Func(hi)) * 7;

	double result = result1 + result2 + result3 + result4;
	result *= step / 90;

	return result;
}

//求積：Romberg公式
double Func_Integral_Romberg
(double lo, double hi, double(*Func)(double), int k = 4)
{//lo:下限；hi:上限；Func:函式；k:等分指數(區間等分成2^k份)
	int size = k + 1;
	double *matrix = new double[size*size];
	for (int i = 0; i < size*size; i++) matrix[i] = 0.0;

	double step = hi - lo;
	matrix[0] = Func_Integral_Trapezoid(lo, hi, Func, 1);
	for (int i = 1; i < size; i++) {
		int n = 1 << (i - 1);
		for (int k = 0; k < n; k++) {
			matrix[i*size + 0] += Func(lo + (k + 0.5)*step);
		}
		matrix[i*size + 0] *= step;
		matrix[i*size + 0] += matrix[(i - 1)*size];
		matrix[i*size + 0] /= 2.0;
		step /= 2.0;
	}

	double temp = 1.0;
	double factor1, factor2;
	for (int j = 1; j < size; j++) {
		temp *= 4.0;
		factor1 = temp / (temp - 1);
		factor2 = 1 / (temp - 1);
		for (int i = j; i < size; i++) {
			matrix[i*size + j] = factor1*matrix[i*size + j - 1]
				- factor2*matrix[(i - 1)*size + j - 1];
		}
	}

	double result = matrix[k*size + k];
	delete[] matrix;

	return result;
}

//測試用的被積函式，0到1積分為PI
double Func_test1(double x)
{
	return 4 / (1 + x*x);
}

int main() {
	cout << "Trapezoid Numerical Integration" << endl;
	cout << setprecision(15) << Func_Integral_Trapezoid(0, 1, Func_test1, 1024) << endl << endl;
	cout << "Simpson Numerical Integration" << endl;
	cout << setprecision(15) << Func_Integral_Simpson(0, 1, Func_test1, 1024) << endl << endl;
	cout << "Cotes Numerical Integration" << endl;
	cout << setprecision(15) << Func_Integral_Cotes(0, 1, Func_test1, 1024) << endl << endl;
	cout << "Romberg Numerical Integraion" << endl;
	cout << setprecision(15) << Func_Integral_Romberg(0, 1, Func_test1, 10) << endl << endl;
	getchar();
	return 0;
}