// Метод Якоби (метод простой итерации) решения СЛАУ
// СЛАУ задается в файле input.txt, ее решение выдается в output.txt
// Допустим, надо решить такую систему:
// 4*x1 -   x2 + 0*x3 = 2
//  -x1 + 4*x2 -   x3 = 6
// 0*x1 -   x2 + 4*x3 = 2
// Тогда файл input.txt должен иметь следующий вид:
//  3
//  4 -1  0  2
// -1  4 -1  6
//  0 -1  4  2

#include <iostream>
#include <fstream>
#include <vector>
using namespace std;

const double epsilon = 0.0001; // требуемая точность

void ReadLinearSystem(string file_name,
                      vector< vector<double> >& linear_system)
{
  linear_system.resize(0);

    ifstream file_linear_system(file_name);

  if (!file_linear_system.is_open()) {
    return;
    }

  int size_linear_system = 0;
  file_linear_system >> size_linear_system;

  vector< double > equation(size_linear_system + 1);
  for (int i = 0; i < size_linear_system; i++) {
    for (int j = 0; j < size_linear_system + 1; j++) {
      file_linear_system >> equation[j];
    }
    linear_system.push_back(equation);
  }

  file_linear_system.close();
}

// Считается норма матрицы размера N,  M
// Предполагается, что матрица непуста
double MatrixNormRow(vector< vector<double> > matrix)
{
  int n = matrix.size();
  int m = matrix[0].size();
  double norm = 0;
  for (int i = 0; i < n; i++) {
    double sum = 0;
    for (int j = 0; j < m - 1; j++) { // Заменить m - 1 на m при копипасте
      sum += fabs(matrix[i][j]);
    }
    if (norm < sum) {
      norm = sum;
    }
  }

  return norm;
}

// Считается норма матрицы размера N,  M
// Предполагается, что матрица непуста
double MatrixNormColumn(vector< vector<double> > matrix)
{
  int n = matrix.size();
  int m = matrix[0].size();
  double norm = 0;
  for (int j = 0; j < m - 1; j++) { // Заменить m - 1 на m при копипасте
    double sum = 0;
    for (int i = 0; i < n; i++) {
      sum += fabs(matrix[i][j]);
    }
    if (norm < sum) {
      norm = sum;
    }
  }

  return norm;
}

double VectorNorm(vector<double> vec)
{
  double norm = 0;
  for (int i = 0; i < vec.size(); i++) {
    double abs_vec_i = fabs(vec[i]);
    if (norm < abs_vec_i) {
      norm = abs_vec_i;
    }
  }
  return norm;
}

bool MissStopCriterion(vector<double> a, vector<double> b)
{
  vector<double> c = a;
  for (int i = 0; i < a.size(); i++) {
    c[i] -= b[i];
  }
  if (VectorNorm(c) < epsilon) {
    return false;
  } else {
    return true;
  }
}

// http://evlm.stuba.sk/~partner2/DBfiles/Numericalanalysis/Linear%20systems%20of%20equations/Jacobi_method_for_systems_of_linear_equations_EN.pdf
// На вход подается система линейных уравнений в виде матрицы размером N, N + 1.
// В части матрицы размером N, N заданы коеффициенты системы. В поледнем
// столбце - свободные члены. Функция возвращает решение системы с помощью
// метода Якоби.
vector<double> JacobiMethod(vector< vector<double> > linear_system)
{
  int linear_system_size = linear_system.size();
  vector<double> solution(linear_system_size);

  // Построение матрицы С и вектора d
  for (int i = 0; i < linear_system_size; i++) {
    int a_ii = linear_system[i][i];
    if (a_ii != 0) {
      linear_system[i][i] = 0;
      linear_system[i][linear_system_size] /= a_ii;
      for (int j = 0; j < linear_system_size; j++) {
        if (i != j) {
          linear_system[i][j] /= -a_ii;
        }
      }
    } else {
      linear_system[i][i] = 1;
      for (int j = 0; j < linear_system_size; j++) {
        if (i != j) {
          linear_system[i][j] *= -1;
        }
      }
    }
  }

  // Проверка сходимости
  bool convergence = true;
  if (MatrixNormRow(linear_system) > 1) {
    if (MatrixNormColumn(linear_system) > 1) {
      convergence = false;
    }
  }

  if (!convergence) {
    return solution;
  }

  // Есть сходимость, запуск итерационного процеса
  vector<double> cur_iteration(linear_system_size);
  for (int i = 0; i < linear_system_size; i++) {
    cur_iteration[i] = linear_system[i][linear_system_size];
  }
  vector<double> prev_iteration(linear_system_size);

  while (MissStopCriterion(cur_iteration, prev_iteration)) {
    prev_iteration = cur_iteration;
    cur_iteration.assign(linear_system_size, 0);
    for (int i = 0; i < linear_system_size; i++) {
      for (int j = 0; j < linear_system_size; j++) {
        cur_iteration[i] += linear_system[i][j] * prev_iteration[j];
      }
      cur_iteration[i] += linear_system[i][linear_system_size];
    }
  }
  solution = cur_iteration;

  return solution;
}

void WriteLinearSystemSolution(string file_name,
                               const vector<double>& linear_system_solution)
{
  ofstream file_linear_system_solution(file_name);
  for (int i = 0; i < linear_system_solution.size(); i++) {
    file_linear_system_solution << "x" << i + 1 << " = "
                                << linear_system_solution[i] << endl;
  }
  file_linear_system_solution.close();
}

int main (void)
{
  vector< vector<double> > linear_system;
    ReadLinearSystem("input.txt", linear_system);
  vector<double> linear_system_solution = JacobiMethod(linear_system);
  WriteLinearSystemSolution("output.txt", linear_system_solution);
    return 0;
}