---

#include <iostream>
#include <vector>
#include <cmath>
#include <ctime>
#include <cstdlib>
#include <iomanip>

using namespace std;

const int MATR_PREC = 1;  // precision when printing matrices
const int PREC = 8;  // precision when printing floating point numbers

/*
 Написать программу, которая решает систему Ax = b указанными в варианте методами
Решить систему с помощью метода Гаусса с выбором ведущего элемента по столбцу и методом квадратного корня (Холецкого), найти определитель матрицы А.
 a[i,i] = 5i
 a[i,j] = 0.1ij
 i != j
 i, j = 1..n
 N = 15
Вектор b задается случайным образом.
Реализовать проверку правильности написанной программы: выводить норму невязки полученного решения, для LU- и QR-разложений выводить значение нормы   или   и т.д. Нормы можно выбирать любые (указать, какие).
Программы должны быть реализованы с максимально возможной экономией памяти и/или количества операций
 */

vector<vector<float> > generate_matr(int n = 15){
    vector<vector<float> > matr(n, vector<float>(n));
    for(int i = 0; i < n; i++){
        matr[i][i] = 5 * (i + 1);
        for(int j = 0; j < n; j++)
            if(i != j)
                matr[i][j] = 0.1 * (i + 1) * (j + 1);
    }
    return matr;
}

void print_matr(vector<vector<float> >& a, vector<float>& b, int max_rows = 10, int max_cols = 10){
    // set fixed notation
    cout << fixed << setprecision(MATR_PREC);
    int n = a.size();
    for(int i = 0; i < min(n, max_rows); i++){
        for(int j = 0; j < min(n, max_cols); j++)
            cout << setw(5) << 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(5) << "..." << " ";
        if(n > max_cols)
            cout << " ...";
        cout << " | ..." << endl;
    }
    cout << endl;
    // back to normal precision
    cout << fixed << setprecision(PREC);
}

// метод для расчета и вывода нормы невязки
// используется максимум-норма
void print_norm(const vector<vector<float> >& a, const vector<float>& b, const vector<float>& x){
    int n = a.size();
    float max_norm = 0;
    for(int i = 0; i < n; i++){
        float sum = 0;
        for(int j = 0; j < n; j++)
            sum += a[i][j] * x[j];
        max_norm = max(max_norm, fabs(sum - b[i]));
//        cout << sum << " - " << b[i] << " = " << fabs(sum - b[i]) << endl;
    }
    cout <<  "Норма невязки(максимум-норма): " << max_norm << endl;
}

// метод для рассчета СЛАУ по методу Гаусса
void gauss(int n, const vector<vector<float> >& const_a, const vector<float>& const_b){
    auto a = const_a;
    auto b = const_b;

    vector<float> x(n);
    // прямой ход
    // k - номер строки
    for(int k = 0; k < n - 1; k++){
        // выбор ведущего элемента по столбцу
        int e_max = k;
        for(int i = k + 1; i < n; i++)
            if(fabs(a[i][k]) > fabs(a[e_max][k]))
                e_max = i;
        // меняем местами k-ую и e_max-ую строки
        swap(a[k], a[e_max]);
        swap(b[k], b[e_max]);
        // i - номер строки, из которой вычитают k-ую строку
        for(int i = k + 1; i < n; i++){
            float m = a[i][k] / a[k][k];
            for(int j = k; j < n; j++)
                a[i][j] -= m * a[k][j];
            b[i] -= m * b[k];
        }
    }
    // обратный ход
    for(int k = n - 1; k >= 0; k--){
        x[k] = b[k];
        for(int i = k + 1; i < n; i++)
            x[k] -= a[k][i] * x[i];
        x[k] /= a[k][k];
    }
    cout << "Метод Гаусса" << endl;
    for(int i = 0; i < n; i++)
        cout << x[i] << " ";
    cout << endl;
    // выводим норму невязки
    print_norm(const_a, const_b, x);
}

void transform_matrix_to_U(int n, vector<vector<float> >& g){
    // найдем матрицу U из прямого хода метода Гаусса
    for(int k = 0; k < n - 1; k++){
        for(int i = k + 1; i < n; i++){
            float m = g[i][k] / g[k][k];
            for(int j = k; j < n; j++)
                g[i][j] -= m * g[k][j];
        }
    }
}

// метод для рассчета СЛАУ по методу квадратного корня (Холецкого)
// (это возможно только для симметричных положительно определенных матриц)
void cholesky(int n, const vector<vector<float> >& const_a, const vector<float>& const_b){
    // A = G * G^T, где G - нижнетреугольная матрица
    auto g = const_a;
    // найдем матрицу U из прямого хода метода Гаусса
    transform_matrix_to_U(n, g);
    // домножаем элементы на корни диагональных
    for(int i = 0; i < n; i++){
        float sqrt_elem = sqrt(g[i][i]);
        for(int j = 0; j < n; j++)
            g[i][j] /= sqrt_elem;
    }
    // транспонируем матрицу
    for(int i = 0; i < n; i++)
        for(int j = i + 1; j < n; j++)
            swap(g[i][j], g[j][i]);
    // решаем систему G * y = b методом Гаусса
    vector<float> y(n);
    y[0] = const_b[0] / g[0][0];
    for(int i = 1; i < n; i++){
        float sum = 0;
        for(int j = 0; j < i; j++)
            sum += g[i][j] * y[j];
        y[i] = (const_b[i] - sum) / g[i][i];
    }
    // решаем систему G^T * x = y методом Гаусса
    vector<float> x(n);
    x[n-1] = y[n-1] / g[n-1][n-1];
    for(int i = n - 2; i >= 0; i--){
        float sum = 0;
        for(int j = i + 1; j < n; j++)
            sum += g[j][i] * x[j];
        x[i] = (y[i] - sum) / g[i][i];
    }
    cout << "Метод квадратного корня (Холецкого)" << endl;
    for(int i = 0; i < n; i++)
        cout << x[i] << " ";
    cout << endl;
    // выводим норму невязки
    print_norm(const_a, const_b, x);
}

// метод для рассчета определителя матрицы
// воспользуемся преобразованием матрицы к верхнетреугольному виду
float det(int n, const vector<vector<float> >& const_a){
    auto matr = const_a;
    transform_matrix_to_U(n, matr);
    float det = 1;
    for(int i = 0; i < n; i++)
        det *= matr[i][i];
    return det;
}


int main(){
    srand(time(nullptr)); // NOLINT(cert-msc51-cpp)
    cout << fixed << setprecision(PREC);

    // исходные векторы, не изменяются, только копируются
    int n = 15;
    vector<vector<float> > a = generate_matr(n);
    vector<float> b(n);
    for(int i = 0; i < n; i++)
        b[i] = rand() % 10000; // NOLINT(cert-msc50-cpp,cppcoreguidelines-narrowing-conversions)

    cout << "Исходная СЛАУ(обрезана до 10 строк/столбцов):" << endl;
    print_matr(a, b);

    // метод Гаусса
    auto start_time = clock();
    gauss(n, a, b);
    cout << "Итоговое время: " << (double)(clock() - start_time) / CLOCKS_PER_SEC << endl << endl;

    // метод квадратного корня (Холецкого)
    start_time = clock();
    cholesky(n, a, b);
    cout << "Итоговое время: " << (double)(clock() - start_time) / CLOCKS_PER_SEC << endl << endl;

    // считаем определитель матрицы
    // должно выйти 10911135640586240917504 (вроде)
    float perfect_det = 10911135640586240917504.;
    float r_det = det(n, a);
    cout << "Определитель матрицы: " << r_det << endl;
    cout << "Точное значение определителя: " << perfect_det << endl;
    cout << "Абсолютная погрешность: " << fabs(r_det - perfect_det) << endl;
    cout << "Относительная погрешность: " << fabs(r_det - perfect_det) / perfect_det * 100 << "%" << endl;

    return 0;
}