Untitled

std::map<rational, vectorspace> math::find_eigenspaces(matrix mat, std::vector<rational> eigenvalues)
{
    std::map<rational, vectorspace> eigenspaces{};
    rational coef{};
    if (mat.get_size()[0] != mat.get_size()[1])
        throw std::exception("Ошибка при поиске собственных подпространств: матрица должна быть квадратной");
    for (int i = 0; i < static_cast<int>(eigenvalues.size()); i++)
    {
        matrix diag_mat = matrix::create_diag_mat(mat.get_size()[0], eigenvalues[i]);
        matrix mat_to_solve = (mat - diag_mat).concat_columns(matrix(mat.get_size()[0], 1));
        // Преобразуем матрицу A-lE = 0 Гауссом и удаляем пустые строки
        matrix gaussian_elim_mat = gaussian_elimination(mat_to_solve, coef).del_empty_rows();

        std::vector<point> basis_vectors{};
        int m = gaussian_elim_mat.get_size()[0];
        int n = gaussian_elim_mat.get_size()[1];

        // Находим индексы базисных и свободных переменных
        std::vector<int> basis_var_indexes = {};
        std::vector<int> free_var_indexes = {};
        for (int j = 0; j < m; j++)
            for (int k = j; k < n - 1; k++)
                if (gaussian_elim_mat[j][k] != 0)
                {
                    basis_var_indexes.push_back(k);
                    break;
                }
        for (int j = 0; j < mat.get_size()[0]; j++)
            if (std::find(basis_var_indexes.begin(), basis_var_indexes.end(), j) == basis_var_indexes.end())
                free_var_indexes.push_back(j);

        std::vector<int> row_indexes_arr(m);
        std::iota(std::begin(row_indexes_arr), std::end(row_indexes_arr), 0);
        // Вырезаем матрицу, состоящую только из столбцов базисных переменных
        matrix sliced_mat = gaussian_elim_mat.intersection_slice(row_indexes_arr, basis_var_indexes);

        // В цикле перебираем все свободные переменные, и делаем их по очереди едиинцей
        for (auto& f : free_var_indexes) 
        {
            std::vector<rational> solution(mat.get_size()[0]);
            solution[f] = 1;
            matrix mat_b(m, 1);

            if (m != 0) {
                // Составляем вектор-столбец b, состоящий из столбца с текущей свободной переменной * -1 (перенесли вправо)
                for (int j = 0; j < m; j++)
                    mat_b[j][0] = gaussian_elim_mat[j][f] * -1;

                // Решаем полученную расширенную матрицу
                matrix expand_sliced_mat = sliced_mat.concat_columns(mat_b);
                std::vector<rational> gauss_solution = solve_square_matrix_by_gauss(expand_sliced_mat);

                for (int j = 0; j < static_cast<int>(gauss_solution.size()); j++)
                    solution[basis_var_indexes[j]] = gauss_solution[j];
            }
            basis_vectors.push_back(solution);
        }

        eigenspaces[eigenvalues[i]] = vectorspace(basis_vectors);
    }

    return eigenspaces;
}