#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;

double function(double x)//прописана исходная функция, которую приближаем
{
	return x * x * x;
}

//прописываем функцию для решения трехдиагональной матрицы методом прогонки. По порядку что такое аругменты: размерность матрицы, нижняя диагональ, средняя диагональ, верхняя диагональ, свободный вектор, массив который получаем на выходе(коэфициенты для посторения сплайна)
void progonka(int n, double *A, double *B, double *C, double *F, double *M) {

	double *alpha = new double[n];
	double *beta = new double[n];
	alpha[1] = -C[0] / B[0];
	beta[1] = F[0] / B[0];

	for (int i = 1; i < n; i++)
	{
		alpha[i + 1] = -C[i] / (B[i] + A[i] * alpha[i]);
		beta[i + 1] = (F[i] - A[i] * beta[i]) / (B[i] + A[i] * alpha[i]);
	}
	M[n] = (F[n] - beta[n] * A[n]) / (B[n] + A[n] * alpha[n]);
	for (int i = n - 1; i >= 0; i--)
		M[i] = alpha[i + 1] * M[i + 1] + beta[i + 1];
}

int main()
{
	ofstream fout("Значения.txt");
	ofstream fout1("ЗначенияФ.txt");
	int n, i;
	double x0, z, a, b, dfa, dfb;
	//вводим (количество точек-1), начало отрезка, конец отрезка, произвоную на левом конце, производную на правом конце
	cin >> n;
	cin >> a;
	cin >> b;
	cin >> dfa;
	cin >> dfb;
	double *x = new double[n + 1];
	double *y = new double[n + 1];
	double *h = new double[n + 1];
	double *A = new double[n + 1];
	double *B = new double[n + 1];
	double *C = new double[n + 1];
	double *F = new double[n + 1];
	double *M = new double[n + 1];

	for (i = 0; i <= n; i++)
	{
		//вводим точки для разбиения отрезка(х[i]) и точное значении исходной функции (y[i])
		//массивы заполняются нулями просто так, этого можно не делать
		cin >> x[i];
		cin >> y[i];
		h[i] = 0;
		B[i] = 0;
		C[i] = 0;
		A[i] = 0;
		F[i] = 0;
		M[i] = 0;
	}
	//теперь заполняем массив длин отрезков
	for (i = 1; i <= n; i++)
		h[i] = x[i] - x[i - 1];
	//массив элементов средней диагонали
	for (i = 1; i < n; i++)
		B[i] = (h[i] + h[i + 1]) / 3.0;

	B[0] = h[1] / 3.0;
	B[n] = h[n] / 3.0;

	//массив элементов верхней диагонали
	for (i = 0; i < n; i++)
		C[i] = h[i + 1] / 6.0;
	//массив элементов нижней диагонали
	for (i = 1; i <= n; i++)
		A[i] = h[i] / 6.0;
	//массив свободного вектора
	for (i = 1; i < n; i++)
	{
		F[i] = (y[i + 1] - y[i]) / h[i + 1] - (y[i] - y[i - 1]) / h[i];
	}
	F[0] = (y[1] - y[0]) / h[1] - dfa;
	F[n] = dfb - (y[n] - y[n - 1]) / h[n];
	//выполняем прогонку
	progonka(n, A, B, C, F, M);
	//вводим какую-нибудь точку отрезка, для того чтобы вычислить приближенное значение в ней (если Барахнин попросит)
	cin >> x0;
	i = 1;
	//вычислем значение в этой точке
	while (x[i] < x0) i++;
	z = (M[i - 1] * pow(x[i] - x0, 3.) + M[i] * pow(x0 - x[i - 1], 3.) + (6.0 * y[i - 1] - M[i - 1] * h[i] * h[i])*(x[i] - x0)
		+ (6.0 * y[i] - M[i] * h[i] * h[i])*(x0 - x[i - 1])) / (6.0 * h[i]);
	cout << z << endl;
	//z-значение в этой точке

	//теперь выводим в файл набор значений в 10к равномерно расположенных точек, чтобы потом построить график (можно сделать чуть проще и переписать вывод, тогда графики можно будет строить в экселе и не париться)
	for (x0 = a; x0 <= b; x0 += (b - a) / 10000) {
		i = 1;
		while (x[i] < x0) i++;
		z = (M[i - 1] * pow(x[i] - x0, 3.) + M[i] * pow(x0 - x[i - 1], 3.) + (6. * y[i - 1] - M[i - 1] * h[i] * h[i])*(x[i] - x0)
			+ (6. * y[i] - M[i] * h[i] * h[i])*(x0 - x[i - 1])) / (6. * h[i]);
		fout << z << endl;
	}

	//так же в файл выводим точные знаения функции в этих же точках
	for (x0 = a; x0 <= b; x0 += (b - a) / 10000) {
		z = function(x0);
		fout1 << z << endl;
	}

	fout.close();
	fout1.close();
	system("pause");
	return 0;
}