#include<iostream>
#include<iomanip>
#include<conio.h>
#include<cmath>
#include<ctime>
#define zloty = 0.618033988749

double interpolacja_liniowa();
double interpolacja_kwadratowa();
double interpolacja_Lagrangea();
double Cubic_Spline();
double Wmaciez_CS(int size, double * v, double * u, double * w, int zmiana);//funkcja do licznia wyznacznika maciezy do Cubic_Spline() zmiana(0)-głowny wyznacznik,zmiana(1)-kolumny podmienione na v(i)
void mciez_troj(int size, double ** maciez);

//rozwiazywanie rownan
double metoda_bisekcji(double(*F)(double));
double metoda_zgadywana(double(*F)(double));//iteraqcyjna
double zloty_podzial_min(double(*F)(double));
double zloty_podzial(double(*F)(double));
double metoda_falsi(double(*F)(double));
double metoda_siecznych(double(*F)(double));
double metoda_newtona_rapshona(double(*F)(double), double(*FP)(double));
double metoda_newtona_rapshona_2_stopnia(double(*F)(double), double(*FP)(double), double(*FPP)(double));
double metoda_halleya(double(*F)(double), double(*FP)(double), double(*FPP)(double));

//funkcje matematyczne
double funkcja(double x);
double funkcja1(double x);
double funkcja2(double x);
double funkcja3(double x);
double funkcja4(double x);
double funkcja4P(double x);
double funkcja4PP(double x);
double funkcja5(double x);
double funkcja6(double x);

//calki
double motoda_trapezow(double(*F)(double));
double metoda_simpsona(double(*F)(double));
double metoda_monte_carlo(double(*F)(double));//losowa

int main(int argc, char* argv[]) {
	do {
		//std::cout << "interpolacja_liniow: " << interpolacja_liniowa();
		//std::cout << "ninterpolacja_kwadratowa: " << interpolacja_kwadratowa();
		//std::cout << "ninterpolacja_Lagrangea: " << interpolacja_Lagrangea();
		//std::cout << "nCubic_Spline: " << Cubic_Spline();
		//std::cout << "nnZloty Podzial min: " << zloty_podzial_min(funkcja2);
		//std::cout << "nnZloty Podzial: " << zloty_podzial(funkcja);
		//std::cout << "nnmetoda bisekcji: " << metoda_bisekcji(funkcja2);
		//std::cout << "nnmetoda zgadywana: " << metoda_zgadywana(funkcja1);
		//std::cout << "nnmetoda_falsi: " << metoda_falsi(funkcja4);
		//std::cout << "nnmetoda_siecznych: " << metoda_siecznych(funkcja4);
		//std::cout << "nnmetoda_newtona_rapshona: " << metoda_newtona_rapshona(funkcja4, funkcja4P);
		//std::cout << "nnmetoda_newtona_rapshona_2_stopnia: " << metoda_newtona_rapshona_2_stopnia(funkcja4, funkcja4P,funkcja4PP);
		//std::cout << "nnmetoda_halleya: " << metoda_halleya(funkcja4, funkcja4P, funkcja4PP);
		//std::cout << "nnmotoda_trapezow: " << motoda_trapezow(funkcja6);
		//std::cout << "nnmotoda_trapezow: " << metoda_simpsona(funkcja6);
		//std::cout << "nnmetoda_monte_carlo: " << metoda_monte_carlo(funkcja6);
		..std::cout << "nn";
		std::cout << "end";
		if (getch() == 27)break;
		std::cout << "nn";
	} while (1);
	return 0;
}

double interpolacja_liniowa()
{

	double x[3];
	double y[3];

	std::cout << "podaj punkt :n x: ";
	std::cin >> x[2];
	//std::cout << "y: ";
	//std::cin >> y[2];

	std::cout << "podaj punkt 0 :n x0: ";
	std::cin >> x[0];
	std::cout << "y0: ";
	std::cin >> y[0];

	std::cout << "podaj punkt 1 :n x1: ";
	std::cin >> x[1];
	std::cout << "y1: ";
	std::cin >> y[1];


	return y[2] = y[0] + ((y[1] - y[0]) / (x[1] - x[0]))*(x[2] - x[0]);
}

double interpolacja_kwadratowa()
{
	double x[4];
	double y[4];

	std::cout << "podaj punkt :n x: ";
	std::cin >> x[3];
	//std::cout << "y: ";
	//std::cin >> y[3];

	std::cout << "podaj punkt 0 :n x0: ";
	std::cin >> x[0];
	std::cout << "y0: ";
	std::cin >> y[0];

	std::cout << "podaj punkt 1 :n x1: ";
	std::cin >> x[1];
	std::cout << "y1: ";
	std::cin >> y[1];

	std::cout << "podaj punkt 2 :n x2: ";
	std::cin >> x[2];
	std::cout << "y2: ";
	std::cin >> y[2];

	return y[3] = ((x[3] - x[1]) / (x[0] - x[1]))*((x[3] - x[2]) / (x[0] - x[2]))*y[0] + ((x[3] - x[0]) / (x[1] - x[0]))*((x[3] - x[2]) / (x[1] - x[2]))*y[1] + ((x[3] - x[0]) / (x[2] - x[0]))*((x[3] - x[1]) / (x[2] - x[1]))*y[2];
}

double interpolacja_Lagrangea()
{
	int _size;
	std::cout << "ile :";
	std::cin >> _size;

	double *x = new double[_size + 1];

	double *y = new double[_size];
	double *s = new double[_size];
	std::cout << "podaj x do obliczen";
	std::cin >> x[_size];

	for (int i = 0; i < _size; i++) {
		std::cout << "podaj x" << i << ": ";
		std::cin >> x[i];
		std::cout << "podaj y" << i << ": ";
		std::cin >> y[i];
	}

	for (int i = 0; i < _size; i++)
	{
		s[i] = y[i];
		for (int j = 0; j < _size; j++)
		{
			if (x[i] != x[j])
			{
				s[i] *= (x[_size] - x[j]);
			}
		}
		for (int j = 0; j < _size; j++)
		{

			if (x[i] != x[j])
			{

				s[i] /= (x[i] - x[j]);

			}
		}
	}

	double suma = 0;
	for (int i = 0; i < _size; i++)
	{
		suma += s[i];

	}

	return suma;
	delete[] x,y,s;
}

double Cubic_Spline()
{
	int _size = 0;
	std::cout << "ile :";
	std::cin >> _size;
	double mainx = 0;
	double mainy = 0;
	double W;
	std::cout << "dla jakiego x :";
	std::cin >> mainx;
	double *x = new double[_size];
	double *y = new double[_size];
	double *a = new double[_size];
	double *b = new double[_size];
	double *c = new double[_size];
	double *d = new double[_size];
	double *h = new double[_size];
	double *v = new double[_size];
	double *u = new double[_size];
	double *w = new double[_size];



	c[0] = 0;
	c[_size - 1] = 0;
	for (int i = 0;i < _size;i++)
	{
		std::cout << "podaj x" << i << " :";
		std::cin >> x[i];
		std::cout << "podaj y" << i << " :";
		std::cin >> y[i];
	}
	for (int i = 0; i < _size; i++)
	{
		a[i] = y[i];
		//std::cout << "a[" << i << "]" << a[i] << " ";
	}
	//std::cout << "nn";
	for (int i = 0; i < _size - 1; i++)
	{
		h[i] = x[i + 1] - x[i];
		//std::cout << "h[" << i << "]" << h[i] << " ";
	}
	//std::cout << "nn";
	for (int i = 1; i < _size - 1; i++)
	{
		u[i] = (h[i - 1] / (h[i - 1] + h[i]));
		//std::cout << "u[" << i << "]" << u[i] <<" ";
	}
	std::cout << "nn";
	for (int i = 1; i < _size - 1; i++)
	{
		w[i] = (h[i] / (h[i - 1] + h[i]));
		//std::cout << "w[" << i << "]" << w[i] << " ";
	}
	std::cout << "nn";
	for (int i = 1; i < _size - 1; ++i) 
	{
		v[i] = (((y[i + 1] - y[i] / h[i]) - (y[i] - y[i - 1] / h[i - 1])) / (h[i - 1] + h[i]));
		//std::cout << "v[" << i << "]" << v[i] << " ";
	}

	W = Wmaciez_CS(_size,v,u,w,0);
	//std::cout << "nn wyznacznik = " << W << "nn";

	for (unsigned int i = 1; i < _size - 1; ++i) 
	{
		c[i] = (Wmaciez_CS(_size, v, u, w, i) / W);
		//std::cout << "c[" << i << "]" << c[i] << " ";
	}
	//std::cout << "nn";

	//for (unsigned int i = 1; i < _size - 1; ++i)
	//{
	//	std::cout << "w[" << i << "]" << Wmaciez_CS(_size, v, u, w, i) << " ";
	//}
	//std::cout << "nn";

	for (int i = 0; i < _size - 1; ++i) 
	{
		d[i] = ((c[i + 1] - c[i]) / 3 * h[i]);
		//std::cout << "d[" << i << "]" << d[i] << " ";
	}
	//std::cout << "nn";

	for (unsigned int i = 0; i < _size; i++) 
	{
		b[i] = (1 / h[i]) * (y[i + 1] - y[i] - c[i] * (h[i] * h[i]) - d[i] * (h[i] * h[i] * h[i]));

		if (i == _size - 1) 
		{
			b[i] = b[i - 1] + 2 * (c[i - 1] * h[i - 1]) + 3 * (d[i - 1] * h[i - 1]);
		}
		//std::cout << "b[" << i << "]" << b[i] << " ";
	}
	std::cout << "nn";
	//x mniejsze
	if (mainx < x[0])
	{
		mainy = a[0] + b[0] * (mainx - x[0]);
	}
	//x wieksze
	if (mainx > x[_size-1])
	{
		mainy = a[_size - 1] + b[_size-2] * (mainx - x[_size - 1]);
	}
	//x w przedziale
	for (int i = 0; i < _size - 1; i++)
	{
		if (x[i] <= mainx && mainx <= x[i + 1])
		{
			mainy = a[i] + b[i] * (mainx - x[i]) + c[i] * ((mainx - x[i])*(mainx - x[i])) + d[i] * ((mainx - x[i])*(mainx - x[i])*(mainx - x[i]));
		}
	}

	delete[] x,y,a,b,c,d,w,v,u,h;

	return mainy;
	
}

double Wmaciez_CS(int size, double * v, double * u, double * w, int zmiana)
{
	double **maciez = new double *[size];
	//wyznacznik
	double W = 1.0;		
	size -= 2;
	for (unsigned int i = 0; i < size; i++) 
	{
		maciez[i] = new double[size];
	}
	for (unsigned int i = 0; i < size; i++)
	{
		for (unsigned int j = 0; j < size; j++) {

			if (i == j) maciez[i][j] = 2;
			else if (i + 1 == j) maciez[i][j] = u[j];
			else if (i == j + 1) maciez[i][j] = w[j + 1];
			else maciez[i][j] = 0;
		}
		if (zmiana != 0) 
		{
			for (unsigned int i = 0; i < size; i++)
			{
				maciez[i][zmiana - 1] = 3 * v[i + 1];
			}
		}
	}
	for (int w = 0; w < size; w++) {

		mciez_troj(size, maciez);//licze maciez trójkątną
	}
	for (int i = 0; i < size; i++)
	{
		W *= maciez[i][i];
	}
	
	for (unsigned int i = 0; i < size; i++)
	{
		delete[] maciez[i];
	}
	delete[] maciez;

	return W;
}

void mciez_troj(int size, double ** maciez)
{
	for (int i = 0; i < size - 1; i++)
	{
		for (int j = i + 1; j < size; j++)
		{
			for (int q = size - 1; q >= i; q--)
			{
				maciez[j][q] = maciez[j][q] - maciez[j][i] * maciez[i][q] / maciez[i][i];
				}
		}
	}

}

double metoda_bisekcji(double(*F)(double))
{
	double a,b;
	double E = 0.000001;
	double i = 0;
	std::cout << "podaj a :";
	std::cin >> a;
	std::cout << "podaj b :";
	std::cin >> b;
	std::cout<< F(a) << " " << F(b) << "nn";
	std::cout << F(a) * F(b) << "nn";
	if (F(a) * F(b) > 0 /*|| (F(a) * F(b) != -INFINITY)*/)
	{
		while (i<100 && F(a) * F(b) > 0)
		{
			b -= b / 100;
			//std::cout << "F(b) = " << F(b) << "n";
			i++;
		}
		if (F(a) * F(b) > 0)
		{
			std::cout << "nie posiada miejsca zerowego!!!!nFunkcja nie jest ciałga na podanym przedzialenZwracam zeronn";
			return 0.;
		}
		i = 0;
	}

	while(b-a>E && i<60)
	{
		if (F((a + b) / 2) * F(a) < 0)
			b = (a + b) / 2;
		else
		{
			a = (a + b) / 2;
		};
		i++;
	}

	std::cout << "ilosc = " << i <<"nn" ;
	return (a+b)/2;
}

double metoda_zgadywana(double(*F)(double))
{
	double a, b;
	double E = 0.000001;
	double i = 0;
	std::cout << "podaj a :";
	std::cin >> a;
	std::cout << "podaj b :";
	std::cin >> b;
	if (F(a) * F(b) > 0 /*|| (F(a) * F(b) != -INFINITY)*/)
	{
		while (i<100 && F(a) * F(b) > 0)
		{
			b -= b / 100;
			//std::cout << "F(b) = " << F(b) << "n";
			i++;
		}
		if (F(a) * F(b) > 0)
		{
			std::cout << "nie posiada miejsca zerowego!!!!nFunkcja nie jest ciałga na podanym przedzialenZwracam zeronn";
			return 0.;
		}
		i = 0;
	}
	while (F(a) * F(b) < 0)
	{
		a += (b - a) / 1000;
	}

	return a;
}

double zloty_podzial_min(double(*F)(double))
{
	double a,b;

	std::cout << "podaj a :";
	std::cin >> a;
	std::cout << "podaj b :";
	std::cin >> b;
	double k = (sqrt(5.) - 1.) / 2.;
	double E = 0.00000001;
	int  i = 0;
	if (F(a) * F(b) > 0 /*|| (F(a) * F(b) != -INFINITY)*/)
	{
		while (i<100 && F(a) * F(b) > 0)
		{
			b -= b / 100;
			//std::cout << "F(b) = " << F(b) << "n";
			i++;
		}
		if (F(a) * F(b) > 0)
		{
			std::cout << "nie posiada miejsca zerowego!!!!nFunkcja nie jest ciałga na podanym przedzialenZwracam zeronn";
			return 0.;
		}
		i = 0;
	}

	while (b - a>E && i<60)
	{
		if (F(a + (b - a)*k) * F(a) < 0)
			b = a + (b - a)*k;
		else
		{
			a = a + (b - a)*k;
		};
		i++;
	}

	std::cout << "ilosc = " << i << "nn";
	return (a + b) / 2;
}

double zloty_podzial(double(*F)(double))
{
	double a, b;

	std::cout << "podaj a :";
	std::cin >> a;
	std::cout << "podaj b :";
	std::cin >> b;
	if (F(a) * F(b) > 0)
	{
		std::cout << "nie posiada miejsca zerowego!!!!nFunkcja nie jest ciałga na podanym przedzialenZwracam zeronn";
		return 0.;
	}
	else
	{
		double k = (sqrt(5.) - 1.) / 2.;
		double xL = b - k * (b - a);
		double xR = a + k * (b - a);
		double E = 0.0000001;
		std::cout << "F(xl) = " << F(xL) << " f(xr) = " << F(xR) << "n";
		std::cout << "xl = " << xL << " xr = " << xR << "n";
		std::cout << "a = " << a << " b = " << b << "n";
		//std::cout << "xl = " << F(xL) << " xr = " << F(xR) << "n";
		while (b - a > E)//dokladność
		{
			std::cout << "F(xl) = " << F(xL) << " f(xr) = " << F(xR) << "n";
			std::cout << "xl = " << xL << " xr = " << xR << "n";
			std::cout << "a = " << a << " b = " << b << "n";
			if (F(xL) < F(xR))
			{
				// wybierz przedział [a, xR]
				b = xR;
				xR = xL;
				xL = b - k * (b - a);
			}
			else
			{
				// wybierz przedział [xL, b]
				a = xL;
				xL = xR;
				xR = a + k * (b - a);
			}
			std::cout << "xl = " << xL << " xr = " << xR << "n";
		}
		return F((a + b) / 2.);
	}
}

double metoda_falsi(double(*F)(double))
{
	double x0, f0, fp;
	double a, b;
	double E = 0.0000001;

	std::cout << "podaj a :";
	std::cin >> a;
	std::cout << "podaj b :";
	std::cin >> b;
	do
	{
		x0 = (a*F(b) - b*F(a)) / (F(b) - F(a));
		f0 = F(x0);
		fp = F(a);
		if (f0 * fp < 0) b = x0;
		else a = x0;
	} while (fabs(f0) > E);

	return x0;
}

double metoda_siecznych(double(*F)(double))
{
	double x0, f0, fp, fk;
	double a, b;
	double E = 0.0000001;

	std::cout << "podaj a :";
	std::cin >> a;
	std::cout << "podaj b :";
	std::cin >> b;

	fp = F(a);
	fk = F(b);

	do
	{
		x0 = (a*fk - b*fp) / (fk - fp);
		f0 = F(x0);
		a = b; 
		fp = fk;
		b = x0; 
		fk = f0;
	} while (fabs(f0) > E);

	return x0;
}

double metoda_newtona_rapshona(double(*F)(double), double(*FP)(double))
{
	double a, b, x0, x1;
	double E = 0.0000001;
	std::cout << "podaj a :";
	std::cin >> a;
	std::cout << "podaj b :";
	std::cin >> b;

	if (F(a)*F(b) >= 0)
		std::cout << "zly przedzial!!!";
	else
	{
		std::cout << "podaj miejsce startowe: ";
		std::cin >> x0;
		while (fabs(F(x0)) > E )
		{
			x1 = x0 - F(x0) / FP(x0);
			x0 = x1;
		}
	}


	return x0;
}

double metoda_newtona_rapshona_2_stopnia(double(*F)(double), double(*FP)(double), double(*FPP)(double))
{
	double a, b, x0, x1;
	double E = 0.0000001;
	std::cout << "podaj a :";
	std::cin >> a;
	std::cout << "podaj b :";
	std::cin >> b;

	if (F(a)*F(b) >= 0)
		std::cout << "zly przedzial!!!";
	else
	{
		std::cout << "podaj miejsce startowe: ";
		std::cin >> x0;
		while (fabs(F(x0)) > E)
		{
			x1 = x0 - ((FP(x0) - sqrt(pow(FP(x0), 2) - 2 * F(x0)*FPP(x0))) / FPP(x0));
			x0 = x1;
		}
	}

	return x0;
}

double metoda_halleya(double(*F)(double), double(*FP)(double), double(*FPP)(double))
{
	double a, b, x0, x1;
	double E = 0.0000001;
	std::cout << "podaj a :";
	std::cin >> a;
	std::cout << "podaj b :";
	std::cin >> b;


	if (F(a)*F(b) >= 0)
		std::cout << "zly przedzial!!!";
	else
	{
		std::cout << "podaj miejsce startowe: ";
		std::cin >> x0;
		std::cout << x0 << " ";
		while (fabs(F(x0)) > E)
		{
			x1 = x0 - (2 * F(x0)*FP(x0)) / (2 * pow(FP(x0), 2) - F(F(x0))*FPP(x0));
			x0 = x1;
			std::cout << x0 << " ";
		}
	}

	return x0;
}


double funkcja(double x)
{
	if (x == 0)return -1.7e-308;
	else return x - (4 / pow(x, 2));
}

double funkcja1(double x)
{
	if (x == 2)return -1.7e-308;
	else return (pow(x,2)-3)/(x-2);
}

double funkcja2(double x)
{
	return pow(x, 2) - 4 + sin(x);
}

double funkcja3(double x)
{
	return x*x;
}

double funkcja4(double x)
{
	return pow(x,3)+pow(x,2)-3*x-3;
}

double funkcja4P(double x)
{
	return 3*pow(x,2)+2*x-3;
}

double funkcja4PP(double x)
{
	return 6 * x, +2;
}

double funkcja5(double x)
{
	return pow(x, 2) - x + 4;
}

double funkcja6(double x)
{
	return pow(x, 2) + 2*x ;
}

double motoda_trapezow(double(*F)(double))
{
	double a, b, ile,k;
	std::cout << "podaj poczatek przedziału: ";
	std::cin >> a;
	std::cout << "podaj koniec przedziału: ";
	std::cin >> b;
	std::cout << "podaj ile: ";
	std::cin >> ile;
	double suma = 0.;
	k =  (b - a) / ile;

	for (int i = 1; i < ile; i++)
	{
		suma += F(a + i * k);
	}
	suma = (suma + (F(a) + F(b)) / 2) * k;
	return suma;
}

double metoda_simpsona(double(*F)(double))
{
	double a, b, ile, k, x,s;
	std::cout << "podaj poczatek przedziału: ";
	std::cin >> a;
	std::cout << "podaj koniec przedziału: ";
	std::cin >> b;
	std::cout << "podaj ile: ";
	std::cin >> ile;
	double suma = 0.;
	s = 0.;
	k = (b - a) / ile;

	for (int i = 1; i <= ile; i++)
	{
		x = a + i*k;
		s += F(x - k / 2);
		if (i < ile) suma += F(x);
	}
	suma = k / 6 * (F(a) + F(b) + 2 * s + 4 * s);

	return suma;
}

double metoda_monte_carlo(double(*F)(double))
{
	srand(time(NULL));
	double a, b, ile, k, x, s;
	std::cout << "podaj poczatek przedziału: ";
	std::cin >> a;
	std::cout << "podaj koniec przedziału: ";
	std::cin >> b;
	std::cout << "podaj ile: ";
	std::cin >> ile;
	double suma = 0.;
	s = 0.;
	k = b - a;

	for (int i = 1; i <= ile; i++)
	{
		suma += F(a + ((double)rand() / (double)(RAND_MAX + 1)*k));
	}
	suma = k * suma / ile;

	return suma;
}