---

#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <cstdio>

using namespace std;

const int N = 10;
const double eps = 1e-6;

// region util
template <typename T>
void print_matrix(vector<vector<T>> &A){
    for(int i = 0; i < A.size(); i++){
        for(int j = 0; j < A.size(); j++){
            printf("%10.3lf ", A[i][j]);
        }
        printf("\n");
    }
}

template <typename T>
void print_vector(vector<T> &v){
    for(int i = 0; i < v.size(); i++){
        printf("%10.3lf ", v[i]);
    }
    printf("\n");
}

vector<vector<double>> generate_matrix(int n){
    vector<vector<double>> A(n, vector<double>(n));
    for(int i = 0; i < n; i++){
        for(int j = i; j < n; j++){
            A[i][j] = A[j][i] = rand() % 100;
        }
    }
    return A;
}

// вспомогательные методы
double norm(vector<double> &v){
    double res = 0;
    for(double e: v){
        res += e * e;
    }
    return sqrt(res);
}

vector<double> operator-(vector<double> &a, vector<double> &b){
    vector<double> res(a.size());
    for(int i = 0; i < a.size(); i++){
        res[i] = a[i] - b[i];
    }
    return res;
}

vector<double> operator*(vector<vector<double>> &A, vector<double> &v){
    vector<double> res(v.size());
    for(int i = 0; i < v.size(); i++){
        for(int j = 0; j < v.size(); j++){
            res[i] += A[i][j] * v[j];
        }
    }
    return res;
}

vector<double> operator*(double a, vector<double> &v){
    vector<double> res(v.size());
    for(int i = 0; i < v.size(); i++){
        res[i] = a * v[i];
    }
    return res;
}

vector<vector<double>> operator*(vector<vector<double>> &A, vector<vector<double>> &B){
    vector<vector<double>> res(A.size(), vector<double>(B[0].size()));
    for(int i = 0; i < A.size(); i++){
        for(int j = 0; j < B[0].size(); j++){
            for(int k = 0; k < B.size(); k++){
                res[i][j] += A[i][k] * B[k][j];
            }
        }
    }
    return res;
}

vector<double> operator/(vector<double> &v, double a){
    vector<double> res(v.size());
    for(int i = 0; i < v.size(); i++){
        res[i] = v[i] / a;
    }
    return res;
}

// поиск максимума, возвращает tuple<row, value>
tuple<int, double> max(vector<double> &v){
    int i_max = 0;
    for(int i = 1; i < v.size(); i++){
        if(abs(v[i]) > abs(v[i_max])){
            i_max = i;
        }
    }
    return {i_max, v[i_max]};
}

// поиск максимума, возвращает tuple<row, col, value>
tuple<int, int, double> max(vector<vector<double>> &A){
    int i_max = 0, j_max = 1;
    for(int i = 0; i < A.size(); i++){
        for(int j = i + 1; j < A.size(); j++){
            if(fabs(A[i][j]) > fabs(A[i_max][j_max])){
                i_max = i;
                j_max = j;
            }
        }
    }
    return {i_max, j_max, A[i_max][j_max]};
}

vector<vector<double>> transpose(vector<vector<double>> A){
    for(int i = 0; i < A.size(); i++){
        for(int j = i + 1; j < A[0].size(); j++){
            swap(A[i][j], A[j][i]);
        }
    }
    return A;
}
// endregion

struct Eigen{
    vector<double> values;
    vector<vector<double>> vectors;
};

// метод вращений якоби
Eigen jacobi(vector<vector<double>> A){
    Eigen eigen;
    eigen.vectors = vector<vector<double>>(A.size(), vector<double>(A.size()));
    for(int i = 0; i < A.size(); i++){
        eigen.vectors[i][i] = 1;
    }
    int iter = 0;
    while(true){
        int i_max, j_max;
        double a_max;
        tie(i_max, j_max, a_max) = max(A);
        if(fabs(a_max) < eps){
            break;
        }
        double phi = 0.5 * atan(2 * a_max / (A[i_max][i_max] - A[j_max][j_max]));
        double c = cos(phi), s = sin(phi);
        vector<vector<double>> U(A.size(), vector<double>(A.size()));
        for(int i = 0; i < A.size(); i++){
            U[i][i] = 1;
        }
        U[i_max][i_max] = U[j_max][j_max] = c;
        U[i_max][j_max] = -s;
        U[j_max][i_max] = s;
        auto U_T = transpose(U);
        auto tmp = U_T * A;
        A = tmp * U;
        eigen.vectors = eigen.vectors * U;
        iter++;
    }
    eigen.vectors = transpose(eigen.vectors);
    eigen.values = vector<double>(A.size());
    for(int i = 0; i < A.size(); i++){
        eigen.values[i] = A[i][i];
    }
    printf("iterations: %d\n", iter);
    return eigen;
}

// степенной метод(максимальное по модулю собственное значение)
Eigen power(vector<vector<double>> A){
    Eigen eigen;
    eigen.vectors = vector<vector<double>>(A.size(), vector<double>(A.size()));
    for(int i = 0; i < A.size(); i++){
        eigen.vectors[i][i] = 1;
    }
    int iter = 0;
    double lambda;
    while(true){
        vector<double> v = A * eigen.vectors[0];
        lambda = norm(v);
        v = v / lambda;
        auto tmp = v - eigen.vectors[0];
        auto t1 = A*v;
        auto t2 = lambda*v;
        auto t3 = t1-t2;
        if(norm(t3) < eps){
            break;
        }
        eigen.vectors[0] = v;
        iter++;
    }
    eigen.values = vector<double>(A.size());
    for(int i = 0; i < A.size(); i++){
        eigen.values[i] = lambda;
    }
    printf("iterations: %d\n", iter);
    return eigen;
}


int main(){
    srand(time(0));
    vector<vector<double>> A = generate_matrix(N);
    printf("matrix A:\n");
    print_matrix(A);

    printf("\n------------------------\n\n");

    printf("Jacobi:\n");
    Eigen eigen = jacobi(A);
    printf("eigen values:\n");
    print_vector(eigen.values);
    printf("eigen vectors:\n");
    print_matrix(eigen.vectors);
    printf("Ax - lambda x:\n");
    for(int i = 0; i < eigen.values.size(); i++){
        vector<double> v = A * eigen.vectors[i];
        auto tmp = eigen.values[i] * eigen.vectors[i];
        v = v - tmp;
        printf("%10.8lf\n", norm(v));
    }

    printf("\n------------------------\n\n");

    printf("Power:\n");
    eigen = power(A);
    printf("eigen max value: %f\n", eigen.values[0]);
    printf("eigen max vector:\n");
    print_vector(eigen.vectors[0]);
    printf("Ax - lambda x:\n");
    vector<double> v = A * eigen.vectors[0];
    auto tmp = eigen.values[0] * eigen.vectors[0];
    v = v - tmp;
    printf("%10.8lf\n", norm(v));
    return 0;
}