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

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

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;
}