#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <ostream>
#include <vector>

double f(double x, double y) {
  double T = M_PI * x * y;
  double sin2T = std::pow(std::sin(T), 2);
  double cos2T = std::pow(std::cos(T), 2);
  double y2_x2 = y * y + x * x;
  return -2 * M_PI * std::exp(std::pow(std::sin(T), 2)) *
         (2 * M_PI * y2_x2 * cos2T * sin2T - M_PI * y2_x2 * sin2T +
          M_PI * y2_x2 * cos2T);
}

double u(double x, double y) {
  return std::exp(std::pow(std::sin(M_PI * x * y), 2));
}

std::vector<std::vector<double>> MinimalResiduals(const int nmax = 1000,
                                                  const double _eps = 0.0000001,
                                                  const int _n = 3,
                                                  const int _m = 3) {
  int Nmax = nmax; // максимальное число итераций (не менее 1)
  std::vector<std::vector<double>> v; // сеточная функция ѵ (x, y)
  v.resize(_n + 1);
  for (size_t i = 0; i < _n + 1; ++i) {
    v[i].resize(_m + 1);
    for (size_t j = 0; j < _m + 1; ++j) {
      v[i][j] = 0;
    }
  }
  std::vector<std::vector<double>> r_vec; // вектор невязки
  r_vec.resize(_n + 1);
  for (size_t i = 0; i < _n + 1; ++i) {
    r_vec[i].resize(_m + 1);
    for (size_t j = 0; j < _m + 1; ++j) {
      r_vec[i][j] = 0;
    }
  }
  std::vector<std::vector<double>> ar; // вектор невязки
  ar.resize(_n + 1);
  for (size_t i = 0; i < _n + 1; ++i) {
    ar[i].resize(_m + 1);
    for (size_t j = 0; j < _m + 1; ++j) {
      ar[i][j] = 0;
    }
  }
  double r1 = 0;
  int S = 0;         // счетчик итераций
  double eps = _eps; // заданная точность
  double eps_max = 0; // точность, достигнутая на текущей итерации
  double eps_cur = 0; // для подсчета точности на текущей итерации
  double a2, k2, h2; // ненулевые элементы матрицы (-А)
  const int n = _n, m = _m; // размерность сетки
  double r_max = 0;
  double a = 0;
  double b = 1;
  double c = 0;
  double d = 1; // границы области определения уравнения
  int i, j; // индексы
  double r = 0;
  double v_old; // старое значение преобразуемой компоненты вектора
  double v_new; // новое значение преобразуемой компоненты вектора и
  bool flag = false; // условие остановки
  double h = ((b - a) / n);
  double k = ((d - c) / m);
  double z = 0;
  double z_max = 0;
  double tau_s = 0;
  double num = 0; // tau=num/denominator
  double denominator = 1;
  double x, y;
  h2 = -std::pow((n / (b - a)), 2);
  k2 = -std::pow((m / (d - c)), 2);
  a2 = -2 * (h2 + k2);
  for (int i = 0; i < n + 1; ++i) {
    v[i][0] = u(a + i * h, a);
    v[i][m] = u(a + i * h, b);
  }
  for (int j = 0; j < m + 1; ++j) {
    v[0][j] = u(c, c + j * k);
    v[n][j] = u(d, c + j * k);
  }
  do {
    z = 0;
    z_max = 0;
    eps_max = 0;
    r_max = -20;
    num = 0;
    denominator = 0;
    for (j = 1; j < m; ++j) {
      for (i = n - 1; i > 0; --i) {
        r_vec[i][j] = f(a + i * h, c + j * k) -
                      (h2 * (v[i + 1][j] + v[i - 1][j]) +
                       k2 * (v[i][j + 1] + v[i][j - 1]) + v[i][j] * a2);
      }
    }
    for (j = 1; j < m; ++j) {
      for (i = n - 1; i > 0; --i) {
        ar[i][j] =
            (h2 * (r_vec[i + 1][j] + r_vec[i - 1][j]) +
             k2 * (r_vec[i][j + 1] + r_vec[i][j - 1]) + r_vec[i][j] * a2);
        num += ar[i][j] * r_vec[i][j];
        denominator += ar[i][j] * ar[i][j];
      }
    }
    tau_s = num / denominator;
    for (j = 1; j < m; ++j) {
      for (i = n - 1; i > 0; --i) {
        v_old = v[i][j];
        v_new = v_old + tau_s * r_vec[i][j];
        eps_cur = std::fabs(v_old - v_new);
        if (eps_cur > eps_max) {
          eps_max = eps_cur;
        }
        v[i][j] = v_new;
        r = fabs(f(a + i * h, c + j * k) -
                 (h2 * (v[i + 1][j] + v[i - 1][j]) +
                  k2 * (v[i][j + 1] + v[i][j - 1]) + v_old * a2));
        if (r > r_max)
          r_max = r;
        z = u(a + i * h, c + j * k) - v[i][j];
        if (z > z_max) {
          z_max = z;
          x = i;
          y = j + 1;
        }
      }
    }

    S = S + 1;
    if ((eps_max <= eps) or (S >= Nmax)) {
      flag = true;
    }
  } while (!flag);

  std::ofstream in("./"
                   "param_test.txt"); // подключение текстового файла для записи
  in << "\nСПРАВКА:" << std::endl;
  in << "Разбиений по x: " << n << '\t' << "разбиений по y:" << n << std::endl;
  in << "Максимальное количество шагов = " << Nmax << std::endl;
  in << "Количество выполненных итераций = " << S << std::endl;
  in << "Требуемая точность = " << eps << std::endl;
  in << "Точность на выходе = " << eps_max << std::endl;
  in << "Невязка = " << r_max << std::endl;
  in << "Норма общей погрешности = " << z_max << std::endl;
  in << "Максимальное отклонение в узле x= " << double(x) / double(n)
     << ", y= " << double(y) / double(m) << std::endl;
  in << "в качестве начального приближения использовано нулевое приближение"
     << std::endl;
  return v;
}

int main() {
  int n = 400;
  int m = 400;
  std::vector<std::vector<double>> v(n, std::vector<double>(m, 0));
  v = MinimalResiduals(2e6, 1e-12, n, m);
  n++;
  m++;
  std::ofstream in("./"
                   "data_test.txt"); // подключение текстового файла для записи
  in << "Результат:" << std::endl;
  for (int i = n - 1; i >= 0; --i) {
    for (int j = 0; j < m; ++j)
      in << std::left << std::setw(10) << v[i][j];
    in << "\n";
  }
  std::ofstream in1(
      "./"
      "data_test_u.txt"); // подключение текстового файла для записи
  in1 << "Результат:" << std::endl;
  for (int i = n - 1; i >= 0; --i) {
    for (int j = 0; j < m; ++j)
      in1 << std::left << std::setw(10) << u(1 / double(i), 1 / double(j));
    in1 << "\n";
  }
  std::ofstream in2(
      "./"
      "data_test_u_v.txt"); // подключение текстового файла для записи
  in2 << "Результат:" << std::endl;
  for (int i = n - 1; i >= 0; --i) {
    for (int j = 0; j < m; ++j)
      in2 << std::left << std::setw(15)
          << std::fabs(u(double(i) / n, double(j) / m) - v[i][j]);

    in2 << "\n";
  }
  return 0;
}