Last active
April 26, 2021 06:59
-
-
Save StefanoBelli/0e9406bae1173119a24a0ae2e740ff82 to your computer and use it in GitHub Desktop.
More than simple math (math utils for C++)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #ifndef NUMERIC_H | |
| #define NUMERIC_H | |
| #include <initializer_list> | |
| #include <type_traits> | |
| #include <sstream> | |
| #include <string> | |
| #include <exception> | |
| #include <vector> | |
| #include <array> | |
| #include <cmath> | |
| namespace ssynx { | |
| namespace num { | |
| // forwarded declaration | |
| namespace matrix { | |
| template <typename T> | |
| class matrix; | |
| } | |
| namespace vector { | |
| class vector_exception : public std::exception { | |
| public: | |
| const char* what() const noexcept { | |
| return "invalid operation, this exclusively depends on the vector itself and/or related operations"; | |
| } | |
| }; | |
| template<typename T> | |
| class vector { | |
| public: | |
| using vector_type = std::vector<T>; | |
| using value_type = T; | |
| using __InitializerType = std::initializer_list<T>; | |
| explicit vector(__InitializerType list) { | |
| vector_type kev; | |
| for(auto const& elem : list) { | |
| kev.push_back(elem); | |
| } | |
| _Initialization(kev); | |
| } | |
| explicit vector(vector_type const& kev) { | |
| _Initialization(kev); | |
| } | |
| explicit vector(matrix::matrix<T> const& mat) { | |
| if(mat.rows() > 1) | |
| throw vector_exception(); | |
| _Initialization(mat.to_vector()[0]); | |
| } | |
| template<std::size_t FixedArrayElem> | |
| explicit vector(std::array<T, FixedArrayElem> const& arr) { | |
| vector_type kev; | |
| for(auto const& elem : arr) { | |
| kev.push_back(elem); | |
| } | |
| _Initialization(kev); | |
| } | |
| vector operator+(vector const& v) const { | |
| if(v.length() != _length) | |
| throw vector_exception(); | |
| vector_type nvec; | |
| for(std::size_t i = 0; i < _length; ++i) { | |
| nvec.push_back(v[i] + vec[i]); | |
| } | |
| return vector<T>(nvec); | |
| } | |
| vector operator-(vector const& v) const { | |
| if(v.length() != _length) | |
| throw vector_exception(); | |
| vector_type nvec; | |
| for(std::size_t i = 0; i < _length; ++i) { | |
| nvec.push_back(v[i] - vec[i]); | |
| } | |
| return vector<T>(nvec); | |
| } | |
| template<typename ScalarType> | |
| vector operator*(ScalarType a) const { | |
| vector_type nvec; | |
| for(std::size_t i = 0; i < _length; ++i) { | |
| nvec.push_back(a * vec[i]); | |
| } | |
| return vector<T>(nvec); | |
| } | |
| //perform dot product | |
| T operator*(vector v) const { | |
| if(v.length() != _length) | |
| throw vector_exception(); | |
| T scalar = 0; | |
| for(std::size_t i = 0; i < _length; ++i) { | |
| scalar += v[i] * vec[i]; | |
| } | |
| return scalar; | |
| } | |
| bool operator==(vector const& v) const { | |
| if(v.length() != _length) | |
| return false; | |
| for(std::size_t i = 0; i < _length; ++i) { | |
| if(v[i] != vec[i]) | |
| return false; | |
| } | |
| return true; | |
| } | |
| bool operator!=(vector const& v) const { | |
| return !operator==(v); | |
| } | |
| bool operator>(vector const& v) const { | |
| if(v.length() != _length) | |
| return false; | |
| for(std::size_t i = 0; i < _length; ++i) { | |
| if(v[i] >= vec[i]) | |
| return false; | |
| } | |
| return true; | |
| } | |
| bool operator<(vector const& v) const { | |
| if(v.length() != _length) | |
| return false; | |
| for(std::size_t i = 0; i < _length; ++i) { | |
| if(v[i] <= vec[i]) | |
| return false; | |
| } | |
| return true; | |
| } | |
| T operator[](std::size_t idx) const { | |
| return at(idx); | |
| } | |
| T at(std::size_t idx) const { | |
| return vec.at(idx); | |
| } | |
| std::size_t length() const noexcept { | |
| return _length; | |
| } | |
| operator bool() const { | |
| return !nulled(); | |
| } | |
| bool nulled() const { | |
| for(auto const& elem : vec) { | |
| if(elem) | |
| return false; | |
| } | |
| return true; | |
| } | |
| std::string printable() const { | |
| std::stringstream ss; | |
| for(auto const& elem : vec) { | |
| ss << elem << ' '; | |
| } | |
| return ss.str(); | |
| } | |
| vector_type const& to_vector() const { | |
| return vec; | |
| } | |
| private: | |
| void _Initialization(vector_type vek) { | |
| vec = std::move(vek); | |
| _length = vec.size(); | |
| } | |
| vector_type vec; | |
| std::size_t _length; | |
| }; | |
| template <typename T> | |
| std::ostream& operator<<(std::ostream& ost, vector<T> const& v) { | |
| ost << v.printable(); | |
| return ost; | |
| } | |
| template <typename ValueType, typename ScalarType> | |
| vector<ValueType> operator*(ScalarType s, vector<ValueType> const& v) { | |
| return v * s; | |
| } | |
| template <typename T> | |
| vector<T> _basic_fill_vector(std::size_t len, T elem) { | |
| typename vector<T>::vector_type vec; | |
| for(std::size_t i = 0; i < len; ++i) | |
| vec.push_back(elem); | |
| return vector<T>(vec); | |
| } | |
| template <typename T = int> | |
| vector<T> nulled_vector(std::size_t len) { | |
| return _basic_fill_vector(len, 0); | |
| } | |
| template <typename T = int> | |
| vector<T> one_filled_vector(std::size_t len) { | |
| return _basic_fill_vector(len, 1); | |
| } | |
| template <typename T = int> | |
| vector<T> canonical_vector(std::size_t len, std::size_t idx) { | |
| typename vector<T>::vector_type vec; | |
| for(std::size_t i = 0; i < len; ++i) { | |
| vec.push_back(i == idx ? 1 : 0); | |
| } | |
| return vector<T>(vec); | |
| } | |
| using vector_int = vector<int>; | |
| using vector_float = vector<float>; | |
| using vector_long = vector<long>; | |
| using vector_double = vector<double>; | |
| using numeric_vector = vector<double>; | |
| } | |
| namespace matrix { | |
| enum class triangular { | |
| NONE, | |
| SUPERIOR, | |
| INFERIOR | |
| }; | |
| class matrix_exception : public std::exception { | |
| public: | |
| const char* what() const noexcept { | |
| return "invalid operation, this exclusively depends on the matrix itself and/or related operations"; | |
| } | |
| }; | |
| namespace misc { | |
| template <typename T> | |
| typename matrix<T>::matrix_type get_new_matrix_support(std::size_t rows, std::size_t cols) { | |
| typename matrix<T>::matrix_type new_matrix_support; | |
| new_matrix_support.resize(rows, typename matrix<T>::vector_type()); | |
| for(auto& col : new_matrix_support) | |
| col.resize(cols); | |
| return new_matrix_support; | |
| } | |
| } | |
| template <typename T> | |
| class matrix { | |
| public: | |
| using __InitializerType = std::initializer_list<std::initializer_list<T>>; | |
| using vector_type = std::vector<T>; | |
| using matrix_type = std::vector<std::vector<T>>; | |
| using value_type = T; | |
| explicit matrix(__InitializerType init_list) { | |
| matrix_type mt; | |
| for(auto const& elem : init_list) | |
| mt.emplace_back(elem); | |
| _Initialization(mt); | |
| } | |
| //may throw if multiple rows have more or less elem | |
| explicit matrix(matrix_type const& vecs) { | |
| _Initialization(vecs); | |
| } | |
| explicit matrix(vector::vector<T> const& numvec) { | |
| _Initialization(matrix_type { numvec.to_vector() }); | |
| } | |
| template<std::size_t RowLength, std::size_t ColLength> | |
| explicit matrix(std::array<std::array<T, ColLength>, RowLength> const& mat) { | |
| auto support = misc::get_new_matrix_support<T>(RowLength, ColLength); | |
| for(std::size_t i = 0; i < RowLength; ++i) { | |
| for(std::size_t j = 0; j < ColLength; ++j) { | |
| support[i][j] = mat[i][j]; | |
| } | |
| } | |
| _Initialization(support); | |
| } | |
| T const& at(std::size_t i, std::size_t j) const { | |
| return _mat.at(i).at(j); | |
| } | |
| std::size_t const& rows() const { | |
| return _rows; | |
| } | |
| std::size_t const& cols() const { | |
| return _cols; | |
| } | |
| matrix col_at(std::size_t j) const { | |
| matrix_type column; | |
| for(auto const& row : _mat) { | |
| column.push_back(vector_type { row[j] } ); | |
| } | |
| return matrix(column); | |
| } | |
| matrix row_at(std::size_t i) const { | |
| return matrix(matrix_type({ _mat[i] })); | |
| } | |
| #define matrix_comparable(m1, action) \ | |
| std::size_t mat1_rows = m1.rows(); \ | |
| std::size_t mat1_cols = m1.cols(); \ | |
| if(mat1_rows != _rows || mat1_cols != _cols) \ | |
| action; \ | |
| //may throw if columns and rows of matrices are not the same | |
| matrix operator+(matrix const& mat1) { | |
| matrix_comparable(mat1, raise_failure()); | |
| matrix_type new_matrix; | |
| for(std::size_t row_idx = 0; row_idx < mat1_rows; ++row_idx) { | |
| new_matrix.emplace_back(vector_type(mat1_cols)); | |
| for(std::size_t col_idx = 0; col_idx < mat1_cols; ++col_idx) { | |
| new_matrix[row_idx][col_idx] = mat1.at(row_idx, col_idx) + _mat[row_idx][col_idx]; | |
| } | |
| } | |
| return matrix(new_matrix); | |
| } | |
| //same as above | |
| matrix operator-(matrix const& mat1) { | |
| matrix_comparable(mat1, raise_failure()); | |
| matrix_type new_matrix; | |
| for(std::size_t row_idx = 0; row_idx < mat1_rows; ++row_idx) { | |
| new_matrix.emplace_back(vector_type(mat1_cols)); | |
| for(std::size_t col_idx = 0; col_idx < mat1_cols; ++col_idx) { | |
| new_matrix[row_idx][col_idx] = mat1.at(row_idx, col_idx) - _mat[row_idx][col_idx]; | |
| } | |
| } | |
| return matrix(new_matrix); | |
| } | |
| //won't throw any matrix_exception | |
| //exception safety is not guranteed, however | |
| template <typename ScalarType> | |
| matrix operator*(ScalarType scalar) const { | |
| matrix_type new_matrix; | |
| for(std::size_t row_idx = 0; row_idx < _rows; ++row_idx) { | |
| new_matrix.emplace_back(vector_type(_cols)); | |
| for(std::size_t col_idx = 0; col_idx < _cols; ++col_idx) { | |
| T good_scalar = std::is_same<T,ScalarType>::value ? scalar : static_cast<T>(scalar); | |
| new_matrix[row_idx][col_idx] = good_scalar * _mat[row_idx][col_idx]; | |
| } | |
| } | |
| return matrix(new_matrix); | |
| } | |
| //may throw if this matrix's cols are different than mat1 rows number (could not perform internal multiplication) | |
| template <typename CompatType> | |
| matrix operator*(matrix<CompatType> const& mat1) { | |
| if(_cols != mat1.rows()) | |
| raise_failure(); | |
| matrix_type resulting_matrix; | |
| resulting_matrix.resize(_rows, vector_type()); | |
| for(std::size_t col_idx = 0; col_idx < mat1.cols(); ++col_idx) { | |
| for(std::size_t row_idx = 0; row_idx < _rows; ++row_idx) { | |
| T component = 0; | |
| for(std::size_t jcol_idx = 0; jcol_idx < _cols; ++jcol_idx) { | |
| component += at(row_idx, jcol_idx) * mat1.at(jcol_idx, col_idx); | |
| } | |
| resulting_matrix[row_idx].push_back(component); | |
| } | |
| } | |
| return matrix(resulting_matrix); | |
| } | |
| //eij(k) | |
| template<typename ScalarType> | |
| matrix elementary_row_mul_and_sum(std::size_t i, std::size_t j, ScalarType k) const { | |
| matrix_type new_matrix = _mat; | |
| vector_type mul_row; | |
| for(auto const& elem : new_matrix[j]) | |
| mul_row.push_back(elem * k); | |
| for(std::size_t idx = 0; idx < _cols; ++idx) | |
| new_matrix[i][idx] += mul_row[idx]; | |
| return matrix(new_matrix); | |
| } | |
| //ei(k) | |
| template<typename ScalarType> | |
| matrix elementary_row_mul_scalar(std::size_t i, ScalarType k) const { | |
| matrix_type new_matrix = _mat; | |
| for(auto& elem : new_matrix[i]) | |
| elem *= k; | |
| return matrix(new_matrix); | |
| } | |
| //Pij | |
| matrix elementary_row_swap(std::size_t i, std::size_t j) const { | |
| vector_type sub_vec_row(_mat[j]); | |
| matrix_type new_matrix = _mat; | |
| new_matrix[j] = _mat[i]; | |
| new_matrix[i] = sub_vec_row; | |
| return matrix(new_matrix); | |
| } | |
| operator bool() const { | |
| return !nulled(); | |
| } | |
| #define matrix_compare(cmp, m, m1) \ | |
| for(std::size_t row_idx = 0; row_idx < _rows; ++row_idx) { \ | |
| for(std::size_t col_idx = 0; col_idx < _cols; ++col_idx) { \ | |
| if(m[row_idx][col_idx] cmp m1.at(row_idx,col_idx)) \ | |
| return false; \ | |
| } \ | |
| } | |
| bool operator==(matrix const& mat1) const { | |
| matrix_comparable(mat1, return false); | |
| matrix_compare(!=, _mat, mat1); | |
| return true; | |
| } | |
| bool operator!=(matrix const& mat1) const { | |
| return !(*this == mat1); | |
| } | |
| bool operator<(matrix const& mat1) const { | |
| matrix_comparable(mat1, return false); | |
| matrix_compare(>=, _mat, mat1); | |
| return true; | |
| } | |
| bool operator>(matrix const& mat1) const { | |
| matrix_comparable(mat1, return false); | |
| matrix_compare(<=, _mat, mat1); | |
| return true; | |
| } | |
| #undef matrix_comparable | |
| #undef matrix_compare | |
| bool nulled() const { | |
| for(auto const& row : _mat) { | |
| for(auto const& col : row) { | |
| if(col) | |
| return false; | |
| } | |
| } | |
| return true; | |
| } | |
| triangular triangularity() const { | |
| if(!squared()) | |
| return triangular::NONE; | |
| triangular is = triangular::INFERIOR; | |
| for(std::size_t row_idx = 0; row_idx < _rows; ++row_idx) { | |
| for(std::size_t col_idx = 0; col_idx < _cols; ++col_idx) { | |
| if(row_idx < col_idx && _mat[row_idx][col_idx]) | |
| is = triangular::SUPERIOR; | |
| } | |
| } | |
| if(is == triangular::SUPERIOR) { | |
| for(std::size_t row_idx = 0; row_idx < _rows; ++row_idx) { | |
| for(std::size_t col_idx = 0; col_idx < _cols; ++col_idx) { | |
| if(row_idx > col_idx && _mat[row_idx][col_idx]) | |
| is = triangular::NONE; | |
| } | |
| } | |
| } | |
| return is; | |
| } | |
| bool diagonal() const { | |
| for(std::size_t row_idx = 0; row_idx < _rows; ++row_idx) { | |
| for(std::size_t col_idx = 0; col_idx < _cols; ++col_idx) { | |
| if(row_idx != col_idx && _mat[row_idx][col_idx]) | |
| return false; | |
| } | |
| } | |
| return true; | |
| } | |
| bool squared() const noexcept { | |
| return _rows == _cols; | |
| } | |
| bool identity() const { | |
| if(!squared() || !diagonal()) | |
| return false; | |
| for(std::size_t row_idx = 0; row_idx < _rows; ++row_idx) { | |
| for(std::size_t col_idx = 0; col_idx < _cols; ++col_idx) { | |
| if(row_idx == col_idx && _mat[row_idx][col_idx] != 1) | |
| return false; | |
| } | |
| } | |
| return true; | |
| } | |
| // *** | |
| std::string printable() const { | |
| std::stringstream printable_matrix; | |
| for (auto const &row : _mat) { | |
| for (auto const &cols : row) { | |
| printable_matrix << cols << ' '; | |
| } | |
| printable_matrix << '\n'; | |
| } | |
| auto s = printable_matrix.str(); | |
| return s.erase(s.size() - 1); | |
| } | |
| matrix_type const& to_vector() const { | |
| return _mat; | |
| } | |
| template <typename _T> | |
| friend long det(matrix<_T> const&); | |
| template <typename _T> | |
| friend long rank(matrix<_T> const&); | |
| template <typename _T> | |
| friend matrix<_T> traspose(matrix<_T> const&); | |
| template <typename _T> | |
| friend matrix<_T> gauss_eliminate(matrix<_T> const&, std::size_t*); | |
| template <typename _T> | |
| friend matrix<_T> align(matrix<_T> const&, matrix<_T> const&); | |
| template <typename _T> | |
| friend matrix<_T> inverse(matrix<_T> const&); | |
| private: | |
| matrix_type _mat; | |
| std::size_t _rows; | |
| std::size_t _cols; | |
| void raise_failure() { | |
| _mat = matrix_type(); | |
| throw matrix_exception(); | |
| } | |
| void _Initialization(matrix_type const& t) { | |
| _mat = t; | |
| _rows = _mat.size(); | |
| _cols = _mat.at(0).size(); | |
| for(auto const& row : _mat) { | |
| if(row.size() != _cols) | |
| raise_failure(); | |
| } | |
| } | |
| }; | |
| namespace misc { | |
| template <typename T> | |
| std::size_t find_column_pivot(matrix<T> const& target, std::size_t start_from) { | |
| for(std::size_t idx = start_from; idx < target.rows(); ++idx) | |
| if(target.at(idx,0)) | |
| return idx; | |
| return 0; | |
| } | |
| //a + by = 0 (b) | |
| //b = -a/y | |
| inline double solve_linear_equation_aplusby(double a, double y) { | |
| if(y == 0) | |
| return -a; | |
| return -a/y; | |
| } | |
| template <typename T> | |
| std::size_t check_if_nulled_col(matrix<T> const& target, std::size_t start_from_row, std::size_t col) { | |
| for(std::size_t row = start_from_row; row < target.rows(); ++row) { | |
| if(target.at(row,col)) | |
| return false; | |
| } | |
| return true; | |
| } | |
| //Ax = 1 | |
| //x = 1/A | |
| inline double scalar_inverse(double a) { | |
| return 1/a; | |
| } | |
| } | |
| template <typename T> | |
| matrix<T> gauss_eliminate(matrix<T> const& mat_, std::size_t* swap_count = nullptr) { | |
| if(mat_._rows == 1) | |
| return mat_; | |
| if(!mat_) | |
| return mat_; | |
| matrix<T> mat = mat_; | |
| std::size_t row_swaps = 0; | |
| for(std::size_t rowidx = 0, colidx = 0; rowidx < mat._rows - 1; ++colidx, ++rowidx) { | |
| auto col = mat.col_at(colidx); | |
| if(misc::check_if_nulled_col(mat, rowidx, colidx)) | |
| continue; | |
| auto pivot = misc::find_column_pivot(col, rowidx); | |
| if(pivot != rowidx) { | |
| mat = mat.elementary_row_swap(pivot, rowidx); | |
| pivot = rowidx; | |
| ++row_swaps; | |
| } | |
| for(std::size_t row = rowidx + 1; row < mat._rows; ++row) { | |
| double k = misc::solve_linear_equation_aplusby(mat.at(row,colidx), mat.at(pivot, colidx)); | |
| if(k == 0) | |
| continue; | |
| mat = mat.elementary_row_mul_and_sum(row, pivot, k); | |
| } | |
| } | |
| if(swap_count) | |
| *swap_count = row_swaps; | |
| return mat; | |
| } | |
| template<typename T> | |
| std::ostream& operator<<(std::ostream& prev_ost, matrix<T> const& m) { | |
| prev_ost << m.printable(); | |
| return prev_ost; | |
| } | |
| template<typename ValueType, typename ScalarType> | |
| matrix<ValueType> operator*(ScalarType s, matrix<ValueType> const& m) { | |
| return m * s; | |
| } | |
| template <typename T = int> | |
| matrix<T> identity_matrix(std::size_t order) { | |
| typename matrix<T>::matrix_type mat; | |
| mat.resize(order, typename matrix<T>::vector_type()); | |
| for(std::size_t row = 0; row < order; ++row) { | |
| mat[row].resize(order); | |
| for(std::size_t col = 0; col < order; ++col) { | |
| if(row == col) | |
| mat[row][col] = 1; | |
| } | |
| } | |
| return matrix<T>(mat); | |
| } | |
| template <typename T = int> | |
| matrix<T> nulled_matrix(std::size_t m, std::size_t n) { | |
| return matrix<T>(misc::get_new_matrix_support<T>(m, n)); | |
| } | |
| template <typename T = int> | |
| matrix<T> one_filled_matrix(std::size_t m, std::size_t n) { | |
| auto support = misc::get_new_matrix_support<T>(m, n); | |
| for(auto& row : support) { | |
| for(auto& elem : row) { | |
| elem = 1; | |
| } | |
| } | |
| return matrix<T>(support); | |
| } | |
| template <typename T = int> | |
| matrix<T> canonical_matrix(std::size_t m, std::size_t n, std::size_t ith, std::size_t jth) { | |
| auto support = misc::get_new_matrix_support<T>(m, n); | |
| support[ith][jth] = 1; | |
| return matrix<T>(support); | |
| } | |
| //may throw if "squared_mat" is not squared | |
| template <typename T> | |
| long det(matrix<T> const& squared_mat) { | |
| if(!squared_mat.squared()) | |
| throw matrix_exception(); | |
| if(squared_mat._rows == 1) | |
| return squared_mat.at(0,0); | |
| if(squared_mat._rows == 2) | |
| return (squared_mat.at(0,0) * squared_mat.at(1,1)) - (squared_mat.at(1,0) * squared_mat.at(0,1)); | |
| long determinant = -1; | |
| std::size_t rowswaps = 0; | |
| auto _eliminated = gauss_eliminate(squared_mat, &rowswaps); | |
| determinant = std::pow(determinant, rowswaps); | |
| for(std::size_t k = 0; k < squared_mat._rows; ++k) | |
| determinant *= _eliminated.at(k,k); | |
| return determinant; | |
| } | |
| template <typename T> | |
| long rank(matrix<T> const& mat) { | |
| if(mat._rows == 0) | |
| return 0; | |
| if(mat._rows == 1) { | |
| if(mat.row_at(0)) | |
| return 1; | |
| else | |
| return 0; | |
| } | |
| long counter = 0; | |
| auto eliminated = gauss_eliminate(mat ); | |
| for(std::size_t row = 0; row < eliminated._rows; ++row) { | |
| if(eliminated.row_at(row)) | |
| counter++; | |
| } | |
| return counter; | |
| } | |
| //mat col to row | |
| template <typename T> | |
| matrix<T> traspose(matrix<T> const& mat) { | |
| typename matrix<T>::matrix_type new_matrix; | |
| new_matrix.resize(mat._cols, typename matrix<T>::vector_type()); | |
| for(std::size_t rowidx = 0; rowidx < mat._cols; ++rowidx) { | |
| new_matrix[rowidx].resize(mat._rows); | |
| for(std::size_t colidx = 0; colidx < mat._rows; ++colidx) | |
| new_matrix[rowidx][colidx] = mat.at(colidx, rowidx); | |
| } | |
| return matrix<T>(new_matrix); | |
| } | |
| template <typename T> | |
| matrix<T> align(matrix<T> const& mat, matrix<T> const& mat1) { | |
| if(mat._rows != mat1._rows) | |
| throw matrix_exception(); | |
| std::size_t total_cols = mat._cols + mat1._cols; | |
| typename matrix<T>::matrix_type new_matrix; | |
| new_matrix.resize(mat._rows, typename matrix<T>::vector_type()); | |
| for(std::size_t rowidx = 0; rowidx < mat._rows; ++rowidx) { | |
| new_matrix[rowidx].resize(total_cols); | |
| for(std::size_t colidx = 0; colidx < mat._cols; ++colidx) | |
| new_matrix[rowidx][colidx] = mat.at(rowidx, colidx); | |
| for(std::size_t colidx = 0; colidx < mat1._cols; ++colidx) | |
| new_matrix[rowidx][colidx+mat1._cols] = mat1.at(rowidx, colidx); | |
| } | |
| return matrix<T>(new_matrix); | |
| } | |
| template<typename T> | |
| matrix<T> inverse(matrix<T> const& mat) { | |
| if(!mat.squared()) | |
| throw matrix_exception(); | |
| if(det(mat) == 0) | |
| return matrix<T>({{}}); //B matrix such that AB = BA = I does not exist | |
| auto intermediary = align(mat, identity_matrix<T>(mat._rows)); | |
| intermediary = gauss_eliminate(intermediary); | |
| for(std::size_t i = 0; i < mat._rows; ++i) { | |
| double k = misc::scalar_inverse(intermediary.at(i, i)); | |
| intermediary = intermediary.elementary_row_mul_scalar(i, k); | |
| } | |
| for (std::size_t j = mat._cols - 1; j >= 0; --j) { | |
| if (j == -1) break; | |
| for (std::size_t i = mat._rows - 1; i >= 0; --i) { | |
| if (i == -1) break; | |
| if (j > i) { | |
| auto k = misc::solve_linear_equation_aplusby(intermediary.at(i, j), 1.0); | |
| intermediary = intermediary.elementary_row_mul_and_sum(i, j, k); | |
| } | |
| } | |
| } | |
| auto support = misc::get_new_matrix_support<T>(mat._rows, mat._cols); | |
| for(std::size_t i = 0; i < mat._rows; ++i) { | |
| for(std::size_t j = 0; j < mat._cols; ++j) { | |
| support[i][j] = intermediary.at(i, mat._cols + j ); | |
| } | |
| } | |
| return matrix<T>(support); | |
| } | |
| std::string to_string(triangular tri) { | |
| if(tri == triangular::NONE) | |
| return "none"; | |
| else if(tri == triangular::SUPERIOR) | |
| return "superior"; | |
| return "inferior"; | |
| } | |
| using matrix_int = matrix<int>; | |
| using matrix_long = matrix<long>; | |
| using matrix_float = matrix<float>; | |
| using matrix_double = matrix<double>; | |
| using numeric_matrix = matrix_double; | |
| } | |
| // vector one-row summatory | |
| template<typename T> | |
| T sum(std::vector<T> const& vek, std::size_t start = 0, std::size_t end = std::string::npos) { | |
| std::size_t ending = end; | |
| if(end == std::string::npos) | |
| ending = vek.size(); | |
| if(start >= end) | |
| return 0; | |
| T __sum = 0; | |
| for(std::size_t i = start; i < ending; ++i) | |
| __sum += vek[i]; | |
| return __sum; | |
| } | |
| template<typename T, std::size_t ArraySize> | |
| T sum(std::array<T, ArraySize> const& arr, std::size_t start = 0, std::size_t end = std::string::npos) { | |
| std::size_t ending = end; | |
| if(end == std::string::npos) | |
| ending = ArraySize; | |
| if(start >= end) | |
| return 0; | |
| T __sum = 0; | |
| for(std::size_t i = start; i < ending; ++i) | |
| __sum += arr[i]; | |
| return __sum; | |
| } | |
| template<typename T> | |
| T sum(vector::vector<T> const& vek, std::size_t start = 0, std::size_t end = std::string::npos) { | |
| return sum(vek.to_vector(), start, end); | |
| } | |
| template<typename T> | |
| T prod(std::vector<T> const& vek, std::size_t start = 0, std::size_t end = std::string::npos) { | |
| std::size_t ending = end; | |
| if(end == std::string::npos) | |
| ending = vek.size(); | |
| if(start >= end) | |
| return 0; | |
| T __prod = 1; | |
| for(std::size_t i = start; i < ending; ++i) | |
| __prod *= vek[i]; | |
| return __prod; | |
| } | |
| template<typename T, std::size_t ArraySize> | |
| T prod(std::array<T, ArraySize> const& arr, std::size_t start = 0, std::size_t end = std::string::npos) { | |
| std::size_t ending = end; | |
| if(end == std::string::npos) | |
| ending = ArraySize; | |
| if(start >= end) | |
| return 0; | |
| T __prod = 1; | |
| for(std::size_t i = start; i < ending; ++i) | |
| __prod *= arr[i]; | |
| return __prod; | |
| } | |
| template<typename T> | |
| T prod(vector::vector<T> const& vek, std::size_t start = 0, std::size_t end = std::string::npos) { | |
| return prod(vek.to_vector(), start, end); | |
| } | |
| template<typename T> | |
| T sum(std::vector<std::vector<T>> const& mat, std::size_t is = 0, std::size_t ie = std::string::npos, | |
| std::size_t js = 0, std::size_t je = std::string::npos) { | |
| std::size_t re = ie; | |
| if(ie == std::string::npos) | |
| re = mat.size(); | |
| T __sum = 0; | |
| for(std::size_t i = is; i < re; ++i) { | |
| std::size_t colsize = je; | |
| if(colsize == std::string::npos) | |
| colsize = mat[i].size(); | |
| for(std::size_t j = js; j < colsize; ++j) { | |
| __sum += mat[i][j]; | |
| } | |
| } | |
| return __sum; | |
| } | |
| template<typename T, std::size_t RowSize, std::size_t ColSize> | |
| T sum(std::array<std::array<T, ColSize>, RowSize> const& mat, std::size_t is = 0, std::size_t ie = std::string::npos, | |
| std::size_t js = 0, std::size_t je = std::string::npos) { | |
| std::size_t re = ie; | |
| if(ie == std::string::npos) | |
| re = RowSize; | |
| std::size_t colsize = je; | |
| if(colsize == std::string::npos) | |
| colsize = ColSize; | |
| T __sum = 0; | |
| for(std::size_t i = is; i < re; ++i) { | |
| for(std::size_t j = js; j < colsize; ++j) { | |
| __sum += mat[i][j]; | |
| } | |
| } | |
| return __sum; | |
| } | |
| template<typename T> | |
| T sum(matrix::matrix<T> const& mat, std::size_t is = 0, std::size_t ie = std::string::npos, | |
| std::size_t js = 0, std::size_t je = std::string::npos) { | |
| return sum(mat.to_vector(), is, ie, js, je); | |
| } | |
| template<typename T> | |
| T prod(std::vector<std::vector<T>> const& mat, std::size_t is = 0, std::size_t ie = std::string::npos, | |
| std::size_t js = 0, std::size_t je = std::string::npos) { | |
| std::size_t re = ie; | |
| if(ie == std::string::npos) | |
| re = mat.size(); | |
| T __prod = 1; | |
| for(std::size_t i = is; i < re; ++i) { | |
| std::size_t colsize = je; | |
| if(colsize == std::string::npos) | |
| colsize = mat[i].size(); | |
| for(std::size_t j = js; j < colsize; ++j) { | |
| __prod *= mat[i][j]; | |
| } | |
| } | |
| return __prod; | |
| } | |
| template<typename T, std::size_t RowSize, std::size_t ColSize> | |
| T prod(std::array<std::array<T, ColSize>, RowSize> const& mat, std::size_t is = 0, std::size_t ie = std::string::npos, | |
| std::size_t js = 0, std::size_t je = std::string::npos) { | |
| std::size_t re = ie; | |
| if(ie == std::string::npos) | |
| re = RowSize; | |
| std::size_t colsize = je; | |
| if(colsize == std::string::npos) | |
| colsize = ColSize; | |
| T __prod = 1; | |
| for(std::size_t i = is; i < re; ++i) { | |
| for(std::size_t j = js; j < colsize; ++j) { | |
| __prod *= mat[i][j]; | |
| } | |
| } | |
| return __prod; | |
| } | |
| template<typename T> | |
| T prod(matrix::matrix<T> const& mat, std::size_t is = 0, std::size_t ie = std::string::npos, | |
| std::size_t js = 0, std::size_t je = std::string::npos) { | |
| return prod(mat.to_vector(), is, ie, js, je); | |
| } | |
| template <typename T> | |
| T prod(std::initializer_list<T> init, std::size_t i = 0, std::size_t end = std::string::npos) { | |
| std::vector<T> vek; | |
| for(auto const& elem : init) | |
| vek.push_back(elem); | |
| return prod(vek, i, end); | |
| } | |
| template <typename T> | |
| T prod(std::initializer_list<std::initializer_list<T>> init, std::size_t i = 0, std::size_t end = std::string::npos, | |
| std::size_t j = 0, std::size_t jend = std::string::npos) { | |
| std::vector<std::vector<T>> vek; | |
| for(auto const& elem : init) | |
| vek.emplace_back(elem); | |
| return prod(vek, i, end, j, jend); | |
| } | |
| template <typename T> | |
| T sum(std::initializer_list<T> init, std::size_t i = 0, std::size_t end = std::string::npos) { | |
| std::vector<T> vek; | |
| for(auto const& elem : init) | |
| vek.push_back(elem); | |
| return sum(vek, i, end); | |
| } | |
| template <typename T> | |
| T sum(std::initializer_list<std::initializer_list<T>> init, std::size_t i = 0, std::size_t end = std::string::npos, | |
| std::size_t j = 0, std::size_t jend = std::string::npos) { | |
| std::vector<std::vector<T>> vek; | |
| for(auto const& elem : init) | |
| vek.emplace_back(elem); | |
| return sum(vek, i, end, j, jend); | |
| } | |
| } | |
| } | |
| #endif |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I would suggest to use numeric_matrix for MOST usages