#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double *cholesky(double *A, int n) {
	double *L = (double*)calloc(n * n, sizeof(double));
	if (L == NULL)
		exit(EXIT_FAILURE);

	for (int i = 0; i < n; i++)
	for (int j = 0; j < (i + 1); j++) {
		double s = 0;
		for (int k = 0; k < j; k++)
			s += L[i * n + k] * L[j * n + k];
		L[i * n + j] = (i == j) ?
			sqrt(A[i * n + i] - s) :
			(1.0 / L[j * n + j] * (A[i * n + j] - s));
	}

	return L;
}

int res;

void show_matrix(double *A, int n) {
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++)
			printf("%lf ", A[i*n + j]);
		printf("\n");
	}
}

void changeMatr(double *A, int n, char* file_name) {
	FILE *input = fopen(file_name, "r");
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++)
			fscanf(input, "%f", A+i*n+j);
	}
	fclose(input);
}

int main() {
	int n = 10;

	double* m = new double[n*n];
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			m[i*n + j] = 1.0 / (i + j + 1.0);
		}
	}
	
	double* b = new double[n];
	for (int i = 0; i < n; i++) {
		b[i] = 1;
	}

	show_matrix(m, n);

	double *res = cholesky(m, n);

	printf("---------\n");
	
	show_matrix(res, n);

	double* mT = new double[n*n];
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			mT[i*n + j] = res[j*n + i];
		}
	}

	printf("---------\n");

	show_matrix(mT, n);

	printf("=========\n");
	changeMatr(m, 10, "input1.txt");
	show_matrix(m, n);
	printf("=========\n");
	
	double* y = new double[n];
	double* x = new double[n];

	for (int i = 0; i < n; i++) {
		y[i] = b[i];
		for (int j = 0; j < i; j++) {
			y[i] -= res[i*n + j] * y[j];
		}
		y[i] /= res[i*n + i];
	}

	for (int i = n - 1; i >= 0; i--) {
		x[i] = y[i];
		for (int j = n - 1; j > i; j--)	{
			x[i] -= mT[i*n + j] * x[j];
		}
		x[i] /= mT[i*n + i];
	}

	printf("---------\n");
	for (int i = 0; i < n; i++) {
		printf("[ ");
		for (int j = 0; j < n; j++)
			printf("%lf ", res[i*n + j]);
		printf("] [%lf] = [%lf]\n", y[i], b[i]);
	}
	printf("---------\n");
	for (int i = 0; i < n; i++) {
		printf("[ ");
		for (int j = 0; j < n; j++)
			printf("%lf ", mT[i*n + j]);
		printf("] [%lf] = [%lf]\n", x[i], y[i]);
	}
	printf("---------\n");
	for (int i = 0; i < n; i++) {
		printf("[ ");
		for (int j = 0; j < n; j++)
			printf("%lf ", m[i*n + j]);
		printf("] [%lf] = [%lf]\n", x[i], b[i]);
	}

	free(res);
	free(m);
	free(mT);
	free(b);
	free(x);
	free(y);
	
	//getchar();

	return 0;
}