#include <valarray>
#include <algorithm>
#include <iostream>
#include <iomanip>

using namespace std;

template<class T> class matrix {
  private:
    valarray<T> v;
    size_t m, n;

  public:
    matrix(size_t m_ = 0, size_t n_ = 0, T val_ = T())
      : v(val_, m_ * n_),m(m_),n(n_) {}

    size_t nrows() const { return m; }
    size_t ncolumns() const { return n; }

    matrix<T>& operator=(T x) { v = x; return *this; }

    class matrix_slice : public slice {
      friend class matrix;
      private:
        valarray<T>& v;

        matrix_slice(valarray<T>& v_, size_t start, size_t size, size_t stride)
          : slice(start, size, stride), v(v_) {}

      public:
        T& operator[](size_t i)
          { return v[start() + i*stride()]; }

        matrix_slice& operator=(const matrix_slice&) = delete;
        matrix_slice& operator=(valarray<T> w) {
          v[static_cast<slice>(*this)] = w;
          return *this;
        }
        matrix_slice& operator=(T x) {
          v[static_cast<slice>(*this)] = x;
          return *this;
        }

        operator valarray<T>() {
          return v[static_cast<slice>(*this)];
        }
    };

    matrix_slice operator[](size_t i) { return row(i); }
    matrix_slice row(size_t i)
      { return matrix_slice{v, i*n, n, 1}; }
    matrix_slice column(size_t j)
      { return matrix_slice{v, j, m, n}; }

    matrix_slice diag()
      { return matrix_slice{v, 0, min(m, n), n + 1}; }

    auto begin() { return std::begin(v); }
    auto end() { return std::end(v); }

    friend istream& operator>>(istream& stream, matrix<T>& a) {
      for (T& x: a)
        stream >> x;
      return stream;
    }
};

using Vektor = valarray<double>;
using Matrix = matrix<double>;

bool jacobi(Matrix a, Vektor& ew, Matrix& ev, int maxcyc = 50, double eps = 1e-12) {
  size_t n{a.nrows()};
  double sum_ndq{0};

  for (size_t i = 0; i < n; i++)
    for (size_t j = i + 1; j < n; j++)
      sum_ndq += a[i][j] * a[i][j];
  sum_ndq *= 2;

  double sum_ndq_start{sum_ndq};

  ev = 0;
  ev.diag() = 1;

  bool success = false;
  for (int cyc = 0; cyc < maxcyc; cyc++) {
    for (size_t i = 0; i < n; i++)
      for (size_t j = i + 1; j < n; j++) {
        double c, s, t, tau;
        if (abs(a[i][j]) <= eps) {
          c = 1;
          s = 0;
        } else {
          tau = (a[j][j] - a[i][i]) / (2 * a[i][j]);
          t = (tau >= 0 ? 1 : -1) / (abs(tau) + sqrt(1 + tau * tau));
          c = 1 / sqrt(1 + t * t);
          s = c * t;
        }

        sum_ndq -= 2 * a[i][j] * a[i][j];

        Vektor temp_i{a.column(i)};
        Vektor temp_j{a.column(j)};
        a.column(i) = c * temp_i - s * temp_j;
        a.column(j) = s * temp_i + c * temp_j;

        temp_i = a.row(i);
        temp_j = a.row(j);
        a.row(i) = c * temp_i - s * temp_j;
        a.row(j) = s * temp_i + c * temp_j;

        temp_i = ev.column(i);
        temp_j = ev.column(j);
        ev.column(i) = c * temp_i - s * temp_j;
        ev.column(j) = s * temp_i + c * temp_j;
      }

    if (sum_ndq <= eps * sum_ndq_start) {
      success = true;
      break;
    }
  }

  ew = a.diag();
  return success;
}

int main() {
  size_t n;
  cout << "n: "; cin >> n;

  Matrix a(n, n);
  cout << "Matrix: " << endl;
  cin >> a;

  Matrix ev(n, n);
  Vektor ew(n);
  if (jacobi(a, ew, ev))
    cout << "Jacobi-Verfahren erfolgreich" << endl;
  else
    cout << "Jacobi-Verfahren nicht erfolgreich" << endl;

  cout << "Eigenwerte:" << endl;
  for (double x: ew)
    cout << "  " << setprecision(4) << setw(8) << fixed << x << endl;
  cout << endl;

  cout << "Eigenvektoren:" << endl;
  for (size_t i = 0; i < n; i++) {
    cout << "  ";
    for (size_t j = 0; j < n; j++)
      cout << setprecision(4) << setw(8) << fixed << ev[i][j];
    cout << endl;
  }
}