lab5

#include<bits/stdc++.h>
using namespace std;

struct Matrix {
  vector<vector<double> > A;
  int n,m;
  Matrix(int n, int m = -1) {
    if (m < 0)
      m = n;
    this->n = n;
    this->m = m;
    A.assign(n, vector<double> (m, 0));
  }

  __inline__ vector<double> & operator [] (const size_t i) {
    return A[i];
  }

  void set_E() {
    if(n == m) {
      for(int i=0; i<n; i++) {
        for(int j=0; j<m; j++) {
          A[i][j] = i == j;
        }
      }
    }
  }

  void randomize() {
    for(int i=0; i<n; i++) {
      for(int j=0; j<m; j++) {
        A[i][j] = rand() * 1.0 / RAND_MAX - 0.5;
      }
    }
  }

  Matrix transponse () {
    Matrix ans(m, n);
    for(int i=0; i<n; i++) {
      for(int j=0; j<m; j++) {
        ans[j][i] = A[i][j];
      }
    }
    return ans;
  }

  bool is_symmetric() {
    if(n != m) {
      return false;
    }
    int ans = 0;
    for(int i=0; i<n; i++) {
      for(int j=0; j<m; j++) {
        ans += A[i][j] != A[j][i];
      }
    }
    return ans == 0;
  }

  double mistake() {
    double answer = 0;
    for (int i=0; i<n; i++) {
      for(int j = i + 1; j < m; j++) {
        answer += A[i][j] * A[i][j];
      }
    }
    return sqrt(2 * answer);
  }
};

inline void mult(Matrix & A, Matrix & B, Matrix & C, int n) {
    for(int i=0;i<n;i++) {
        for(int j=0;j<n;j++) {
            C[i][j] = 0;
            for(int k=0;k<n;k++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
}

inline void transp(Matrix & A, Matrix & B, int n) {
    for(int i=0;i<n;i++) {
        for(int j=0;j<n;j++) {
            B[i][j] = A[j][i];
        }
    }
}

Matrix rotating_method(Matrix & A, double eps, size_t & iterations) {
  const double multiplier =  sqrt(2) / 2;
  iterations = 0;
  int n = A.n;
  Matrix RM(n),RMT(n), buff(n), ans(n);
  ans.set_E();
  double max_value = 0;
  int x = 0, y = 0;
  double angle = 0;
  while(A.mistake() > eps) {
    iterations++;
    RM.set_E();
    x = 0, y = 0;
    max_value = 0;
    for(int i=0; i<n; i++) {
      for(int j= i + 1; j<n; j++) {
        if(abs(A[i][j]) > max_value) {
          max_value = abs(A[i][j]);
          x = i, y = j;
        }
      }
    }
    if(A[x][x] == A[y][y]) {
      RM[x][x] = RM[y][y] = RM[y][x] = multiplier;
      RM[x][y] = - multiplier;
    } else {
      angle = 0.5 * atan(2.0 * A[x][y] / (A[x][x] - A[y][y]));
      RM[x][x] = RM[y][y] = cos(angle);
      RM[x][y] = -sin(angle);
      RM[y][x] = sin(angle);
    }
        transp(RM, RMT, n);
    mult(RMT, A, buff, n);
    mult(buff, RM, A, n);
    mult(ans, RM, buff, n);
    ans = buff;
  }
  return ans;
}

int main() {
  ios_base::sync_with_stdio(false);
  cin.tie(nullptr);
  int n;
  cin>>n;
  double eps;
  cin>>eps;
  Matrix A(n);
  A.randomize();
  Matrix At = A.transponse();
  Matrix buff(n);
  mult(A, At, buff, n);
  A = buff;
  double diag = 0;
  for(int i=0; i<n; i++) {
    diag += A[i][i];
  }
  size_t iterations = 0;
  Matrix ans = rotating_method(A, eps, iterations);
  for(int i=0; i<n; i++) {
    cout<<"Eigenvector "<<i + 1<<" : ";
    for(int j=0; j<n; j++) {
      cout<<ans[j][i]<<" ";
    }
    cout<<endl;
  }
  cout<<"Eigenvalues: ";
  double sum_of_eigenvalues = 0;
  for(int i=0; i<n; i++) {
    cout<<A[i][i]<<" ";
    sum_of_eigenvalues += A[i][i];
  }
  cout<<endl;
  cout<<"Total iteration steps : "<<iterations<<endl;
  cout<<"Diagonal sum : "<<diag<<endl;
  cout<<"Eigenvalues sum : "<<sum_of_eigenvalues<<endl;
  return 0;
}