#include <iostream>
#include <stdlib.h>
#include <fstream>
#include <string>
#include <iomanip>
#include <math.h>
double TOLX = 1.0e-5;
double TOL0 = 1.0e-5;
int nmax = 15;
const int SIZE = 4;

double A[SIZE][SIZE] = {	{100.,1.,-2.,3.},
							{4.,300.,-5.,6.},
							{7.,-8.,400.,9.},
							{-10.,11.,-12.,200.} };
double B[SIZE] = { 395.,603.,-415.,-606. };

double estymator(double x2[],double x1[]){
	double max = 0;
	double tab[SIZE];
	for (int i = 0; i < SIZE; i++) {
		tab[i] = x2[i] - x1[i];
		if (fabs(tab[i]) > max)
			max = fabs(tab[i]);
	}
	std::cout << "Estymator: " << max << std::endl;
	return max;
}/*
double residuum(double a[][SIZE],double b[],double x[]){
	double r[SIZE] = {}, max = 0;
	for (int i = 0; i < 4; i++) {
		for (int j = 0; j < 4; j++) {
			r[i] += A[i][j] * x[j];
			std::cout << r[i] << " ";
		}
		r[i] = fabs(r[i] - b[i]);
		if (r[i] > max)
			max = r[i];
	}
	for (int i = 0; i < SIZE; i++)
	{
		if (r[i] >= max)
			max = r[i];
	}
	std::cout << "m ax" << max << "  ";

	return max;
}*/
double reziduum(double x[], double b[], const double a[][SIZE]) {

	double r[SIZE] = { 0.,0.,0.,0. };
	double max;
	for (int i = 0; i < SIZE; i++) {
		for (int j = 0; j < SIZE; j++) {
			r[i] += a[i][j] * x[j];
		}
		std::cout << r[i] << " ";
		r[i] = fabs(r[i] - b[i]);
		std::cout <<"po "<< r[i] << " ";
	}
	max = r[0];
	for (int i = 0; i < SIZE; i++) {
		std::cout << r[i] << " ";
		if (r[i] > max) {
			//std::cout << r[i] << " ";
			max = r[i];
		}
	}
	return max;
}
double jacobi(double arr[][SIZE], double b[SIZE]) {
	double D[SIZE], b2[SIZE], x2[SIZE], x1[SIZE] = { 1.,1.,1.,1. };
	int n = 0;
	//D^-1
	for (int i = 0; i < SIZE; i++)
	{
		D[i] = 1. / arr[i][i];
		arr[i][i] = 0;
		//ok
	}
	//b * D --> wspolczynnik c
	for (int i = 0; i < SIZE; i++)
	{
		b2[i] = D[i] * b[i];
		//ok
	}
	//przeksztalcenie macierzy A
	for (int i = 0; i < SIZE; i++)
	{
		for (int j = 0; j < SIZE; j++)
		{
			if(i!=j)
			arr[i][j] = arr[i][j] * ((-1)*D[i]);
		}
	}
	int num = 0;
	//algorytm
	for (int k = 0; k < nmax; k++)
		{
			for (int i = 0; i < SIZE; i++)
			{
				x2[i] = b2[i];
				for (int j = 0; j < SIZE; j++)
				{
					x2[i] += arr[i][j] * x1[j];
				}
			}
			num++;
			std::cout << "Iteracja " << num << "  ";
			for (int i = 0; i < SIZE; i++)
			{
				std::cout << "x[" << i << "]=" << x2[i] << "  ";
			}
			//residuum(A, B, x2);
			std::cout << "Residuum="<<reziduum(x2, B, A) << std::endl;
			if (estymator(x2, x1) < TOLX && reziduum(x2, B, A) < TOL0)break;
			for (int i = 0; i < SIZE; i++) {
				x1[i] = x2[i];
			}
			std::cout << std::endl<<std::endl;
		}
	return 0;
}
int main() {

	double X[SIZE] = {1.,1.,1.,1.};
	jacobi(A, B);
	getchar();
	return 0;
}