#define _USE_MATH_DEFINES
#include <iostream>
#include <math.h>
#include <fstream>
#include <string>
using namespace std;

double aLeft = 0.0, bRight = M_PI / 4.0;
double alfa = 0.0, beta = 1.0, gamma = 0.0, fi = 0.0, psi = 1.0, theta = -0.5;
double P = 1.0, Q = 0.0, R = 4.0;

double s(double x) {
	return tan(x);
}

double analyticalSolution(double x) {
	return (2.0 * x * cos(2.0 * x) + 2 * sin(2.0 * x) - log(2.0) * sin(2.0 * x) - 2 * log(cos(x)) * sin(2.0 * x)) / 4.0;
}

double getEstimator(double* misscalc, int N) {
	double e = fabs(misscalc[0]);
	for (int i = 1; i < N; i++)
		if (fabs(misscalc[i]) > e) {
			e = fabs(misscalc[i]);
		}
	return e;
}

void thomasAlgorithm(double** A, int N) {
	for (int i = 1; i < N; i++) {
		A[1][i] -= A[2][i - 1] * A[0][i - 1] / A[1][i - 1];
	}
}

void solution(double** A, double* b, double* x, int N) {
	for (int i = 1; i < N; i++) {
		b[i] -= A[2][i - 1] * b[i - 1] / A[1][i - 1];
	}

	x[N - 1] = b[N - 1] / A[1][N - 1];

	for (int i = N - 2; i >= 0; i--) {
		x[i] = (b[i] - A[0][i] * x[i + 1]) / A[1][i];
	}
}

void saveData(double* x, double h, int i, string name) {
	double x_i = aLeft;
	std::fstream file;
	file.open(name, fstream::out);
	file.precision(15);
	for (int j = 0; j < i; j++) {
		file << x_i << " " << x[j] << " " << analyticalSolution(x_i) << " " << endl;
		x_i += h;
	}
	file.close();
}

double conventional(double h, int i) {
	double** A = new double*[3];
	A[0] = new double[i - 1]; //upper
	A[1] = new double[i]; //diagonal
	A[2] = new double[i - 1]; //lower
	double* b = new double[i];
	double* misscalc = new double[i];
	double* x = new double[i];
	double x_i = aLeft, e;

	A[0][0] = alfa / h;
	A[1][0] = beta - alfa / h;
	b[0] = -gamma;

	for (int j = 1; j < i - 1; j++) {
		A[0][j] = (P / h / h) + (Q / 2.0 / h);
		A[1][j] = R - (2.0 * P / h / h);
		A[2][j - 1] = (P / h / h) - (Q / 2.0 / h);
		x_i += h;
		b[j] = -s(x_i);
	}

	A[1][i - 1] = fi / h + psi;
	A[2][i - 2] = -fi / h;
	b[i - 1] = -theta;

	thomasAlgorithm(A, i);
	solution(A, b, x, i);

	x_i = aLeft;

	for (int j = 0; j < i; j++) {
		misscalc[j] = fabs(x[j] - analyticalSolution(x_i));
		x_i += h;
	}
	if (i == 100) {
		saveData(x, h, i, "conventional.txt");
	}
	e = getEstimator(misscalc, i);

	delete[] x;
	delete[] b;
	delete[] misscalc;
	for (int j = 0; j < 3; j++) {
		delete[] A[j];
	}
	delete[] A;

	return e;
}

double Numerow(double h, int i) {
	double** A = new double* [3];
	A[0] = new double[i - 1]; //upper
	A[1] = new double[i]; //diagonal
	A[2] = new double[i - 1]; //lower
	double* b = new double[i];
	double* misscalc = new double[i];
	double* x = new double[i];
	double x_i = aLeft, e;

	A[0][0] = alfa / h;
	A[1][0] = beta - alfa / h;
	b[0] = -gamma;

	for (int j = 1; j < i - 1; j++) {
		A[0][j] = (P / (h * h)) + (R / 12.0);
		A[1][j] = (R * 10.0 / 12.0) - (2.0 * P / (h * h));
		A[2][j - 1] = (P / (h * h)) + (R / 12.0);
		x_i += h;
		b[j] = -(s(x_i - h) / 12.0 + (10.0 / 12.0) * s(x_i) + s(x_i + h) / 12.0);
	}

	A[1][i - 1] = fi / h + psi;
	A[2][i - 2] = -fi / h;
	b[i - 1] = -theta;

	thomasAlgorithm(A, i);
	solution(A, b, x, i);

	x_i = aLeft;
	for (int j = 0; j < i; j++) {
		misscalc[j] = fabs(x[j] - analyticalSolution(x_i));
		x_i += h;
	}

	if (i == 100) {
		saveData(x, h, i, "Numerow.txt");
	}
	e = getEstimator(misscalc, i);

	delete[] x;
	delete[] b;
	delete[] misscalc;
	for (int j = 0; j < 3; j++) {
		delete[] A[j];
	}
	delete[] A;

	return e;
}

int main()
{
	double h;
	fstream misscal;
	misscal.open("misscalculations.txt", ios::out);
	misscal.precision(15);
	h = (bRight - aLeft) / 9.0;
	misscal << log10(h) << " " << log10(conventional(h, 10)) << " " << log10(Numerow(h, 10)) << endl;
	for (int i = 100; i < 100000; i += 100) {
		h = (bRight - aLeft) / (i - 1);
		misscal << log10(h) << " " << log10(conventional(h, i)) << " " << log10(Numerow(h, i)) << endl;
	}
	misscal.close();
}