#include <stdio.h>
#include <stdbool.h>
#include <unistd.h>
#include <stdlib.h>
#include <float.h>
#include <string.h>
#include <math.h>

typedef double matr_elem;
struct matrix {
	unsigned short n;
	unsigned short m;
	matr_elem **elems;
};

/*
 * ==========================================================================
 * Utility functions
 * ==========================================================================
 */

struct matrix *build_matrix
(unsigned short n, unsigned short m, matr_elem elems[])
{
	struct matrix* matr = (struct matrix*)malloc(sizeof(struct matrix));
	matr->n = n;
	matr->m = m;
        matr->elems = (matr_elem**)malloc(m * sizeof(matr_elem*));
	for (unsigned short i = 0; i < m; ++i) {
		matr->elems[i] = (matr_elem*)malloc(n * sizeof(matr_elem));
		for (unsigned short j = 0; j < n; ++j)
			matr->elems[i][j] = elems[i*n + j];
	}
	return matr;
};

void free_matrix (struct matrix *matr)
{
	for (unsigned short i = 0; i < matr->m; ++i)
		free(matr->elems[i]);
	free(matr->elems);
	free(matr);
}

void print_matrix (struct matrix *matr)
{
	if (!matr)
		return;

	printf("(");
	for (unsigned short i = 0; i < matr->m; ++i) {
		if (i != 0)
			printf(" ");
		printf("( ");
		for (unsigned short j = 0; j < matr->n; ++j)
			printf("%.6lf\t", matr->elems[i][j]);
		printf(")");
		if (i == matr->m - 1)
			printf(")");
		printf("\n");
	}
}

/*
 * ==========================================================================
 * Computational functions
 * ==========================================================================
 */

bool diagonal_no_zeroes (struct matrix *a, struct matrix *b)
{
	for (unsigned short i = 0; i < a->m; ++i)
		if (fabs(a->elems[i][i]) < DBL_EPSILON) {
			bool found_nonzero = false;
			for (unsigned short j = 0; (j < a->m) && !found_nonzero; ++i)
				if (fabs(a->elems[j][i]) > DBL_EPSILON) {
					found_nonzero = true;
					for (unsigned short k = 0; k < a->n; ++k)
						a->elems[i][k] += a->elems[j][k];
					b->elems[i][0] += b->elems[j][0];
				}
			if (!found_nonzero)
				return false;
		}
	return true;
}

bool matrix_is_dominant (struct matrix *matr)
{
	bool found_strict = false;
	for (unsigned short i = 0; i < matr->m; ++i) {
		matr_elem sum;
		for (unsigned short j = 0; j < matr->n; ++j)
			sum += abs(matr->elems[i][j]) * (j != i ? 1 : -1);
		if (sum > DBL_EPSILON)
			return false;
		if ((sum > -DBL_EPSILON) && (sum < DBL_EPSILON))
			found_strict = true;
	}
	return found_strict;
}

bool consistent_sle (struct matrix *a)
{
	if (!a)
		return NULL;

	for (int i = 0; i < a->m; ++i)
		for (int j = 0; j < a->n; ++j)
			if ((i != j) &&
				(abs(a->elems[i][j]) >= abs(a->elems[i][i])))
					return false;
	return true;
}

struct matrix *comp_sle (struct matrix *a, struct matrix *b, double eps,
			 unsigned int M)
{
	if ((b->n != 1) || (a->n != a-> m) || (a->m != b->m))
		return NULL;

	matr_elem *x = (matr_elem*)calloc(a->n, sizeof(matr_elem));
	double m;
	unsigned int k = 0;
	do {
		++k;
		m = 0;
		for (unsigned short i = 0; i < a->n; ++i) {
			matr_elem p = x[i];
			double s = 0;
			for (unsigned short j = 0; j < a->n; ++j)
				if (i != j)
					s += a->elems[i][j] * x[j];
			x[i] = (b->elems[i][0] - s) / a->elems[i][i];
			m += (x[i] - p) * (x[i] - p);
		}
	} while ((m > eps*eps) && (k < M));

	struct matrix *matr = build_matrix (1, a->n, x);
	free(x);
	if (k >= M)
		return NULL;
	return matr;
}

/*
 * ==========================================================================
 * Main body
 * ==========================================================================
 */

int main ()
{
	unsigned short m;
	printf("Please input the number of equations: ");
	bool correct = true;
	char **check = 0;
	do {
		char s[1024];
		scanf("%s", s);
		m = strtol(s, check, 10);
		if ((m <= 0) || (**check)) {
			m = 0;
			printf("Неверное m\n");
		}
	} while (!m);
	unsigned short n = m;
	matr_elem *a_elems = (matr_elem*)malloc((n * m) * sizeof(matr_elem));
	matr_elem *b_elems = (matr_elem*)malloc(n * sizeof(matr_elem));
	printf("Input the coefficients and constant terms of each equation"
			"separated by spaces.\n");
	for (unsigned short i = 0; i < m; ++i)
		for (unsigned short j = 0; j <= n; ++j)
		{
			correct = true;
			do {
			char s[1024];
			scanf("%s", s);
			double *p = (j != n ? a_elems+i*m+j : b_elems+i);
			*p = strtof(s, check);
			correct = **check == '\0';
			if (!correct) {
				printf("Неверное число\n");
			}
			} while (!correct);
		}
	struct matrix *a = build_matrix(n, m, a_elems);
	struct matrix *b = build_matrix(1, m, b_elems);
	if (!diagonal_no_zeroes(a, b)) {
		printf("Zeroed column!\n");
		exit(1);
	}
	bool c = consistent_sle(a);
	printf("Matrix is %sdominant", (matrix_is_dominant(a) ? "" : "non-"));
	printf("\n");
	printf("SLE is %sconsistent\n", (c ? "" : "in"));
	free(a_elems);
	free(b_elems);
	double eps;
	printf("Please input the E of the calculation: ");
	do {
		char s[1024];
		scanf("%s", s);
		eps = strtof(s, check);
		correct = **check == '\0';
		if (!correct) {
			printf("Неверное число\n");
		}
	} while (!correct);
	unsigned int M;
	printf("Please input the M of the calculation: ");
	do {
		char s[1024];
		scanf("%s", s);
		char **check = 0;
		M = strtol(s, check, 10);
		if ((M == 0) || (**check)) {
			printf("Неверное M\n");
		}
	} while (!M);
	struct matrix *x = comp_sle(a, b, eps, M);
	free_matrix(a);
	free_matrix(b);
	if (!x) {
		printf("Does not compute.\n");
		exit(1);
	}
	print_matrix(x);
	free_matrix(x);
}