Оптимизировать или сделать код немножко красивее)

[Pein95]

#include <iostream>
#include <conio.h>
#include <cmath>

using namespace std;

double rectangle(double a, double b, double eps);
double trap (double a, double b, double eps);
double simpson (double a, double b, double eps);
double F(double x);

int main()
{
	double square,eps;
	cout << "eps = ";
	cin >> eps;
	square = rectangle(0,3.14/4.0,eps);
	cout << square << '\n';
	square = trap(0,3.14/4.0,eps);
	cout << square << '\n';
	square = simpson(0,3.14/4.0,eps);
	cout << square << '\n';
	getch();
	return 0;
}

double rectangle(double a, double b, double eps)
{


	double func1=0, func2=0;
	int n=1;
	double h;
	double xi;
	
	do
	{
		func1=0;
		func2=0;
		h=(b-a)/n;
		xi=a;
		for (int i=0; i<n; i++)
		{
			func1 +=pow(1+sin(2*xi),2.0);
			xi+=h;
		}
		func1 *=h;
		n*=2;
		xi = a;
		h=(b-a)/n;
		for (int i=0; i<n; i++)
		{
			func2 +=pow(1+sin(2*xi),2.0);
			xi+=h;
		}
		func2 *=h;
		n*=2;
	
	}while(fabs(func2-func1)>eps);
	return func2;
}

double trap(double a, double b, double eps)
{
	
	double func1=0, func2=0;
	int n=1;
	double h;
	double xi;
	double sum=0;

	do
	{
		func1=0;
		func2=0;
		sum=0;
		h = (b-a)/n;
		xi = a;
		func1 = (a+b)/2;
		for (int i=1; i<n; i++)
		{
			xi+=h;
			sum+=pow(1+sin(2*xi),2.0);
		}
		func1 = h*(func1+sum);

		n*=2;
		h = (b-a)/n;
		xi = a;
		sum=0;
		for (int i=1; i<n; i++)
		{
			xi+=h;
			sum+=pow(1+sin(2*xi),2.0);
		}
		func2 = h*(func2+sum);
		
	}while(func2-func1>eps);

	return func2;

}

double simpson (double a, double b, double eps)
{
	
	double h;
	double func1, func2;
	int n=1,i;
	double xi;
	double sum=0;

	do
	{
		func1=0;
		func2=0;
		n*=2;
		h = (b-a)/n;
		i=1;
		xi = a;
		func1 = pow(1+sin(2*xi),2.0);
		do
		{
			xi +=h;
			func1 = func1+4*pow(1+sin(2*xi),2.0);
			i+=2;
			if (!(i>=n))
			{
				xi = xi+h;
				func1 = func1 +2*pow(1+sin(2*xi),2.0);
			}
		}while(i<n);
		xi = b;
		func1 = (func1+pow(1+sin(2*xi),2.0))*(h/3);

		n *=2;
		h = (b-a)/n;
		i=1;
		xi = a;
	
		func2 = pow(1+sin(2*xi),2.0);
		do
		{
			xi +=h;
			func2 = func2+4*pow(1+sin(2*xi),2.0);
			i+=2;
			if (!(i>=n))
			{
				xi = xi+h;
				func2 = func2 +2*pow(1+sin(2*xi),2.0);
			}
		}while(i<n);
		xi = b;
		func2 = (func2+pow(1+sin(2*xi),2.0))*(h/3);


	}while (fabs(func2-func1)>eps);

	return func2;

	

}

double F(double x)
{
	double res = pow(1+sin(2*x),2.0);
	return res;
}

Ну вопшем реализировано три функции, которые щитаю площадь фигуры с точностю eps. Можно ли както оптимизировать ето? Или сделать код более понятным. А то задание простое, а код немножко страшний)