#include <iostream>
#include <cmath>
#include <utility>
#include <vector>

std::vector<std::vector<std::vector<double>>> lu_composition(std::vector<std::vector<double>> matrix_a) {
    std::vector<std::vector<std::vector<double>>> result;

    std::vector<std::vector<double>> matrix_u = std::move(matrix_a);
    std::vector<std::vector<double>> matrix_l = {
            {1, 0, 0, 0},
            {0, 1, 0, 0},
            {0, 0, 1, 0},
            {0, 0, 0, 1}
    };
    for (int root_row = 0; root_row < matrix_u.size(); root_row++)
    {
        double multipler;

        for (int row = root_row + 1; row < matrix_u.size(); row++) {
            for (int col = root_row; col < matrix_u.size(); ++col) {
                if (col == root_row) {
                    multipler = -(matrix_u[row][col] / matrix_u[root_row][col]);
                    matrix_l[row][col] = -multipler;
                }

                matrix_u[row][col] = matrix_u[root_row][col] * multipler + matrix_u[row][col];
            }

//            Displays matrix_u and matrix_l
//
//            for (auto a: matrix_u) {
//                std::cout << "{";
//                for (auto b: a) {
//                    std::cout << b << " ";
//                }
//                std::cout << "}" << std::endl;
//            }
//
//            for (auto a: matrix_l) {
//                std::cout << "{";
//                for (auto b: a) {
//                    std::cout << b << " ";
//                }
//                std::cout << "}" << std::endl;
//            }
//
//            std::cout << std::endl;
        }
    }

    result.push_back(matrix_u);
    result.push_back(matrix_l);

    return result;
};

// std::vector<double>
void forward_subsitution(std::vector<std::vector<double>> matrix_l, std::vector<double> matrix_b) {
    std::vector<double> matrix_y =
            {
                    0,
                    0,
                    0,
                    0
            };

    matrix_y[0] = matrix_b[0];

    for(int i = 1; i < matrix_b.size(); i++) {
        double sum = 0;
        for (int j = 0; j < i - 1; j++) {
            sum =+ matrix_l[i][j] * matrix_y[j];
        }
        matrix_y[i] = (matrix_b[i] - sum) / matrix_l[i][i];
    }


    for(auto variable: matrix_y)
        std::cout << variable << std::endl;
    //return  matrix_y;
}

int main() {
    std::vector<std::vector<double>> matrix_u;
    std::vector<std::vector<double>> matrix_l;
    std::vector<double> matrix_y;
    std::vector<std::vector<double>> matrix_a = {
            {4, -2, 4, 2},
            {3, 1, 4, 2},
            {2, 4, 2, 1},
            {2, -2, 4, 2}
    };

    std::vector<double> matrix_b = {
            8,
            7,
            10,
            2
    };

    std::vector<std::vector<std::vector<double>>> temp = lu_composition(matrix_a);

    matrix_u = std::move(temp[0]);
    matrix_l = std::move(temp[1]);

    forward_subsitution(matrix_l, matrix_b);
}