---

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <iomanip>

using namespace std;

const float eps = 1e-4;
const int kmax = 1000;
const int prec = 5;  // precision when printing floating point numbers
const string sep = "----------------------------------------\n";

// region util
vector<vector<float>> get_matrix(int n){
    vector<vector<float>> matrix(n, vector<float>(n));
    for(int i = 0; i < n; i++){
        matrix[i][i] = 9.f * (float)pow(i+1, 1./3);
        for(int j = 0; j < n; j++)
            if(i != j)
                matrix[i][j] = 0.01f / (float)(sqrt(i + 1) + sqrt(j + 1));
    }
    return matrix;
}

bool check_strict_diagonal_dominance(const vector<vector<float>>& matrix){
    for(int i = 0; i < matrix.size(); i++){
        float sum = 0;
        for(int j = 0; j < matrix.size(); j++)
            if(i != j)
                sum += fabs(matrix[i][j]);
        if(fabs(matrix[i][i]) <= sum)
            return false;
    }
    return true;
}

vector<float> get_solution(const vector<vector<float>>& matrix){
    vector<float> solution(matrix.size());
    for(int i = 0; i < matrix.size(); i++)
        solution[i] = (float)i + 5;
    return solution;
}

vector<float> get_b(const vector<vector<float>>& matrix, const vector<float>& solution){
    vector<float> b(matrix.size());
    for(int i = 0; i < matrix.size(); i++){
        float sum = 0;
        for(int j = 0; j < matrix.size(); j++)
            sum += matrix[i][j] * solution[j];
        b[i] = sum;
    }
    return b;
}

vector<vector<float>> get_modified_matrix(const vector<vector<float>>& matrix){
    auto modified_matrix = matrix;
    for(int i = 1; i < matrix.size(); i++){
        float sum = 0;
        for(int j = 0; j < matrix.size(); j++)
            if(i != j)
                sum += fabs(matrix[i][j]);
        modified_matrix[i][i] = sum;
    }
    return modified_matrix;
}

void print_matr(vector<vector<float> >& a, vector<float>& b, int max_rows = 5, int max_cols = 5, int w = 10){
    // set fixed notation
    cout << fixed << setprecision(prec);
    int n = (int)a.size();
    for(int i = 0; i < min(n, max_rows); i++){
        for(int j = 0; j < min(n, max_cols); j++)
            cout << setw(w) << a[i][j] << " ";
        // if to many columns, print "..."
        if(n > max_cols)
            cout << " ...";
        cout << " | " << b[i] << endl;
    }
    // if to many rows, print "..."
    if(n > max_rows){
        for(int i = 0; i < max_cols; i++)
            cout << setw(w) << "..." << " ";
        if(n > max_cols)
            cout << " ...";
        cout << " | ..." << endl;
    }
    cout << endl;
}

void print_vector(const vector<float>& v, int w = 5){
    cout << fixed << setprecision(prec);
    for(float i: v)
        cout << setw(w) << i << " ";
    cout << endl;
}

float norm(const vector<float>& v){
    float norm = 0;
    for(float i: v)
        norm = max(norm, fabs(i));
    return norm;
}

vector<float> operator-(const vector<float>& a, const vector<float>& b){
    vector<float> c(a.size());
    for(int i = 0; i < a.size(); i++)
        c[i] = a[i] - b[i];
    return c;
}
// endregion

struct solution{
    vector<float> x;
    int k;
};

// якоби
solution jacobi(const vector<vector<float>>& a, const vector<float>& b, [[maybe_unused]] float o = 1){
    int n = (int)a.size();
    vector<float> x(n);
    vector<float> x_prev(n);
    int k = 0;
    do{
        k++;
        x_prev = x;
        for(int i = 0; i < n; i++){
            float sum = 0;
            for(int j = 0; j < n; j++)
                if(i != j)
                    sum += a[i][j] * x_prev[j];
            x[i] = (b[i] - sum) / a[i][i];
        }
    }while(k < kmax && norm(x - x_prev) > eps);
    return {x, k};
}

// гаусс-зейдель
solution gauss_seidel(const vector<vector<float>>& a, const vector<float>& b, [[maybe_unused]] float omega = 1){
    int n = (int)a.size();
    vector<float> x(n);
    vector<float> x_prev(n);
    int k = 0;
    do{
        k++;
        x_prev = x;
        for(int i = 0; i < n; i++){
            float sum = 0;
            for(int j = 0; j < n; j++)
                if(i != j)
                    sum += a[i][j] * (j < i ? x[j] : x_prev[j]);
            x[i] = (b[i] - sum) / a[i][i];
        }
    }while(k < kmax && norm(x - x_prev) > eps);
    return {x, k};
}

// релаксация
solution relaxation(const vector<vector<float>>& a, const vector<float>& b, float omega){
    int n = (int)a.size();
    vector<float> x(n);
    vector<float> x_prev(n);
    int k = 0;
    do{
        k++;
        x_prev = x;
        for(int i = 0; i < n; i++){
            float sum = 0;
            for(int j = 0; j < n; j++)
                if(i != j)
                    sum += a[i][j] * (j < i ? x[j] : x_prev[j]);
            x[i] = (1 - omega) * x_prev[i] + omega * (b[i] - sum) / a[i][i];
        }
    }while(k < kmax && norm(x - x_prev) > eps);
    return {x, k};
}

void test_method(const string& name, solution f(const vector<vector<float>>&, const vector<float>&, float),
                 const vector<vector<float>>& a, const vector<float>& b,
                 const vector<float>& solution, float omega = 1){
    cout << name << ":" << endl;
    auto method_solution = f(a, b, omega);
    print_vector(method_solution.x);
    cout << "Iterations: " << method_solution.k << endl;
    cout << "Is k limit reached: " << (method_solution.k == kmax ? "yes" : "no") << endl;
    float method_norm = norm(solution - method_solution.x);
    float method_relative_error = method_norm / norm(solution);
    cout << "Norm: " << method_norm << endl;
    cout << "Relative error: " << method_relative_error * 100 << "%" << endl;
}

int main(){
    cout << fixed << setprecision(prec);

    int n = 15;
    vector<vector<float>> matrix = get_matrix(n);
    vector<float> solution = get_solution(matrix);
    vector<float> b = get_b(matrix, solution);

    cout << "Matrix:" << endl;
    print_matr(matrix, b);
    cout << "Strict diagonal dominance: ";
    cout << (check_strict_diagonal_dominance(matrix) ? "yes" : "no") << endl;
    cout << "Solution:" << endl;
    print_vector(solution);
    cout << sep;

    test_method("Jacobi", jacobi, matrix, b, solution);
    cout << sep;

    test_method("Gauss-Seidel", gauss_seidel, matrix, b, solution);
    cout << sep;

    test_method("Relaxation (omega = 0.5)", relaxation, matrix, b, solution, 0.5);
    cout << sep;

    test_method("Relaxation (omega = 1.5)", relaxation, matrix, b, solution, 1.5);
    cout << sep;

    auto modified_matrix = get_modified_matrix(matrix);
    auto modified_b = get_b(modified_matrix, solution);
    cout << "Modified matrix:" << endl;
    print_matr(modified_matrix, modified_b);
    cout << "Strict diagonal dominance: ";
    cout << (check_strict_diagonal_dominance(modified_matrix) ? "yes" : "no") << endl;
    cout << "Solution:" << endl;
    print_vector(solution);
    cout << sep;

    test_method("Jacobi", jacobi, modified_matrix, modified_b, solution);
    cout << sep;

    test_method("Gauss-Seidel", gauss_seidel, modified_matrix, modified_b, solution);
    cout << sep;

    test_method("Relaxation (omega = 0.5)", relaxation, modified_matrix, modified_b, solution, 0.5);
    cout << sep;

    test_method("Relaxation (omega = 1.5)", relaxation, modified_matrix, modified_b, solution, 1.5);
    cout << sep;
    return 0;
}