#include "stdafx.h"
#include "stdio.h"
#include "math.h"
#include "conio.h"
#include "iostream.h"
double PI = 3.1415926;

/// Определение функции под интегралом
double function(double x)
{
	return (sin(log(x)))/x; //вот сюда ввёл свою функцию
}

/// интегрирование методом центральных прямоугольников (с шагом)
double middleRectangles(double left, double right, double step)
{
	double I = 0;
	for(double xi = left + step * 0.5; xi < right; xi += step)
	{
		I += function(xi);
	}
	I *= step;
	return I;
}

/// интегрирование методом центральных прямоугольников (по числу итераций)
double middleRectangles(double left, double right, int n)
{
	double I = 0;
	// iin?eoaee oaae
	double step = (right - left) / n;
	for(double xi = left + step * 0.5; xi < right; xi += step)
	{
		I += function(xi);
	}
	I *= step;
	return I;
}

/// интегрирование методом левых прямоугольников (с шагом)
/// y[0] y[n-1]
double leftRectangle(double left, double right, double step)
{
	double I = 0;
	for(double xi = left; xi < right; xi += step)
	{
		I += function(xi);
	}
	I *= step;
	return I;
}
/// интегрирование методом левых прямоугольников(по числу итераций)
/// y[0] y[n-1]
double leftRectangle(double left, double right, int n)
{
	double I = 0;
	double step = (right - left) / (double)n;
	for(double xi = left; xi < right; xi += step)
	{
		I += function(xi);
	}
	I *= step;
	return I;
}

/// интегрирование методом правых прямоугольников
/// y[1] y[n]
double rightRectangle(double left, double right, double step)
{
	double I = 0;
	for(double xi = left + step; xi <= right; xi += step)
	{
		I += function(xi);
	}
	I *= step;
	return I;
}
/// интегрирование методом правых прямоугольников (по числу итераций)
/// y[1] y[n]
double rightRectangle(double left, double right, int n)
{
	double I = 0;
	double step = (right - left) / (double)n;
	for(double xi = left + step; xi <= right; xi += step)
	{
		I += function(xi);
	}
	I *= step;
	return I;
}

/// метод трапеций (по шагу)
/// формула трапеции имеет второй порядок точности (m = 2)
double trapezium(double left, double right, double step)
{
	double I= 0;
	for(double xi = left + step; xi < right; xi += step)
	{
		I += function(xi);
	}
	I = (I + 0.5*(function(left) + function(right))) * step;
	return I;
}

/// метод трапеций (по числу итераций)
/// формула трапеции имеет второй порядок точности (m = 2)
double trapezium(double left, double right, int n)
{
	double I = 0;
	double step = (right - left) / (double)n;
	for(double xi = left + step; xi < right; xi += step)
	{
		I += function(xi);
	}
	I = (I + 0.5*(function(left) + function(right))) * step;
	return I;
}



//---------------------------------------------------------------------------
//метод симпсона итерационный (задается числом итераций)
//порядок точности метода Симпсона = 4
double simpson(double a, double b, int n)
{
	double step = (b - a) / n;
	double I = 0, I1 = 0, I2 = 0;
	// в краевых точках
	double I0 = function(a)+function(b);
	I2 = function(a + step);
	for(int iteration = 2; iteration < n; iteration += 2 )
	{
		// нечетные
		I2 += function(a + (iteration + 1) * step);
		// четные
		I1 += function(a + iteration * step);
	}
	I = I0 + 4*I2 + 2*I1;
	/// домножаем на коэффициент шага
	I *= step/3;

	return I;
}

////// левые прямоугольники с заданной точностью
double leftRectangleEpsilon(double left, double right, double epsilon)
{
	int n = 10;
	double I0 = leftRectangle(left, right, n);
	double I1 = 0;

	/// пока погрешность больше заданной увеличиваем число шагов в два раза
	while( fabs(I0 - I1) > epsilon )
	{
		I1 = I0;
		n *= 2;
		I0 = leftRectangle(left, right, n);
	}
	return I0;
}

/// правые прямоугольники с заданной точностью
double rightRectangleEpsilon(double left, double right, double epsilon)
{
	int n = 10;
	double I0 = rightRectangle(left, right, n);
	double I1 = 0;
	while( fabs(I0 - I1) > epsilon )
	{
		I1 = I0;
		n *= 2;
		I0 = rightRectangle(left, right, n);
	}
	return I0;
}

/// Центральные прямоугольники с заданной точностью
double middleRectangleEpsilon(double left, double right, double epsilon)
{
	int n = 10;
	double I0 = middleRectangles(left, right, n);
	double I1 = 0;
	while( fabs(I0 - I1) > epsilon )
	{
		I1 = I0;
		n *= 2;
		I0 = middleRectangles(left, right, n);
	}
	return I0;
}

/// трапеции с заданной точностью
double trapeziumEpsilon(double left, double right, double epsilon)
{
	int n = 10;
	double I0 = trapezium(left, right, n);
	double I1 = 0;
	while( fabs(I0 - I1) > epsilon )
	{
		I1 = I0;
		n *= 2;
		I0 = trapezium(left, right, n);
	}
	return I0;
}

/// Симпсон с заданной точностью
double simpsonEpsilon(double left, double right, double epsilon)
{
	int n = 10;
	double I0 = simpson(left, right, n);
	double I1 = 0;
	while( fabs(I0 - I1) > epsilon )
	{
		I1 = I0;
		n *= 2;
		I0 = simpson(left, right, n);
	}
	return I0;
}


/// печатаем таблицу
void printTable(double (*Fx)(double, double, int), double left, double right, int n, int step)
{
	cout<<"N          "<<"EPSILON"<<endl;
	for(int i = 2; i <= n; i+=step)
	{
		double I0 = Fx(left, right, i);
		double I1 = Fx(left, right, 2*i);
		double epsilon = fabs(I0-I1);
		cout<<"n = "<<i<<" epsilon = "<< epsilon<<endl;
	}
}

/// оценка погрешности по правилу Рунге
/// |I (h/2) - I(h)|
int main(int argc, char ** argv)
{
	double a, b;
	double epsilon = 0;
	a = 0;
	b = PI / 2;

	/// cin<<a;
	/// cin<<b;

	cout<<"Vuberite"<<endl;
	cout<<"1. Vichislenie integrala"<<endl;
	cout<<"2. Tablici"<<endl;

	int vubor = 2;
	cin>>vubor;

	if(vubor == 2)
	{
		cout<<"Vuberite"<<endl;
		cout<<"1. Levie pryamougolniki"<<endl;
		cout<<"2. Pravie pryamougolniki"<<endl;
		cout<<"3. Centralnie pryamougolniki"<<endl;
		cout<<"4. Simpson"<<endl;
		cout<<"5. Trapecii"<<endl;
		cin>>vubor;

		switch(vubor)
		{
		case 1:
			printTable(leftRectangle, a, b, 30, 2);	
			break;
		case 2:
			printTable(rightRectangle, a, b, 30, 2);
			break;
		case 3:
			printTable(middleRectangles, a, b, 30, 2);
			break;
		case 4:
			printTable(simpson, a, b, 30, 2);
			break;
		case 5:
			printTable(trapezium, a, b, 30, 2);
			break;
		default:
			printTable(trapezium, a, b, 100, 10);
			break;
		}
	}
	else
	{
		printf("Vvedite espilon : ");
		cin>>epsilon;

		double simps = simpsonEpsilon(a, b, epsilon);
		double rectL = leftRectangleEpsilon(a, b, epsilon);
		double rectR = rightRectangleEpsilon(a, b, epsilon);
		double rectM = middleRectangleEpsilon(a, b, epsilon);
		double trap = trapeziumEpsilon(a, b, epsilon);

		printf("Otvet method levih pryamoygolnikov  %10.10f.\n", rectL);
		printf("Otvet method pravih pryamoygolnikov   %10.10f.\n", rectR);
		printf("Otvet method Centralnih pryamoygolnikov   %10.10f.\n", rectM);
		printf("Otvet method simpsona  %10.10f.\n", simps);
		printf("Otvet method trapecii  %10.10f.\n", trap);
	}

	getch();

	return 0;
}