Instantly share code, notes, and snippets.

# StefanoBelli/num.h

Last active April 26, 2021 06:59
Show Gist options
• 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 #include #include #include #include #include #include #include namespace ssynx { namespace num { // forwarded declaration namespace matrix { template 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 class vector { public: using vector_type = std::vector; using value_type = T; using __InitializerType = std::initializer_list; 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 const& mat) { if(mat.rows() > 1) throw vector_exception(); _Initialization(mat.to_vector()[0]); } template explicit vector(std::array 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(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(nvec); } template 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(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 std::ostream& operator<<(std::ostream& ost, vector const& v) { ost << v.printable(); return ost; } template vector operator*(ScalarType s, vector const& v) { return v * s; } template vector _basic_fill_vector(std::size_t len, T elem) { typename vector::vector_type vec; for(std::size_t i = 0; i < len; ++i) vec.push_back(elem); return vector(vec); } template vector nulled_vector(std::size_t len) { return _basic_fill_vector(len, 0); } template vector one_filled_vector(std::size_t len) { return _basic_fill_vector(len, 1); } template vector canonical_vector(std::size_t len, std::size_t idx) { typename vector::vector_type vec; for(std::size_t i = 0; i < len; ++i) { vec.push_back(i == idx ? 1 : 0); } return vector(vec); } using vector_int = vector; using vector_float = vector; using vector_long = vector; using vector_double = vector; using numeric_vector = vector; } 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 matrix::matrix_type get_new_matrix_support(std::size_t rows, std::size_t cols) { typename matrix::matrix_type new_matrix_support; new_matrix_support.resize(rows, typename matrix::vector_type()); for(auto& col : new_matrix_support) col.resize(cols); return new_matrix_support; } } template class matrix { public: using __InitializerType = std::initializer_list>; using vector_type = std::vector; using matrix_type = std::vector>; 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 const& numvec) { _Initialization(matrix_type { numvec.to_vector() }); } template explicit matrix(std::array, RowLength> const& mat) { auto support = misc::get_new_matrix_support(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 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::value ? scalar : static_cast(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 matrix operator*(matrix 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 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 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 friend long det(matrix<_T> const&); template friend long rank(matrix<_T> const&); template friend matrix<_T> traspose(matrix<_T> const&); template friend matrix<_T> gauss_eliminate(matrix<_T> const&, std::size_t*); template friend matrix<_T> align(matrix<_T> const&, matrix<_T> const&); template 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 std::size_t find_column_pivot(matrix 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 std::size_t check_if_nulled_col(matrix 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 matrix gauss_eliminate(matrix const& mat_, std::size_t* swap_count = nullptr) { if(mat_._rows == 1) return mat_; if(!mat_) return mat_; matrix 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 std::ostream& operator<<(std::ostream& prev_ost, matrix const& m) { prev_ost << m.printable(); return prev_ost; } template matrix operator*(ScalarType s, matrix const& m) { return m * s; } template matrix identity_matrix(std::size_t order) { typename matrix::matrix_type mat; mat.resize(order, typename matrix::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(mat); } template matrix nulled_matrix(std::size_t m, std::size_t n) { return matrix(misc::get_new_matrix_support(m, n)); } template matrix one_filled_matrix(std::size_t m, std::size_t n) { auto support = misc::get_new_matrix_support(m, n); for(auto& row : support) { for(auto& elem : row) { elem = 1; } } return matrix(support); } template matrix 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(m, n); support[ith][jth] = 1; return matrix(support); } //may throw if "squared_mat" is not squared template long det(matrix 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 long rank(matrix 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 matrix traspose(matrix const& mat) { typename matrix::matrix_type new_matrix; new_matrix.resize(mat._cols, typename matrix::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(new_matrix); } template matrix align(matrix const& mat, matrix const& mat1) { if(mat._rows != mat1._rows) throw matrix_exception(); std::size_t total_cols = mat._cols + mat1._cols; typename matrix::matrix_type new_matrix; new_matrix.resize(mat._rows, typename matrix::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(new_matrix); } template matrix inverse(matrix const& mat) { if(!mat.squared()) throw matrix_exception(); if(det(mat) == 0) return matrix({{}}); //B matrix such that AB = BA = I does not exist auto intermediary = align(mat, identity_matrix(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(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(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; using matrix_long = matrix; using matrix_float = matrix; using matrix_double = matrix; using numeric_matrix = matrix_double; } // vector one-row summatory template T sum(std::vector 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 T sum(std::array 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 T sum(vector::vector const& vek, std::size_t start = 0, std::size_t end = std::string::npos) { return sum(vek.to_vector(), start, end); } template T prod(std::vector 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 T prod(std::array 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 T prod(vector::vector const& vek, std::size_t start = 0, std::size_t end = std::string::npos) { return prod(vek.to_vector(), start, end); } template T sum(std::vector> 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 T sum(std::array, 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 T sum(matrix::matrix 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 T prod(std::vector> 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 T prod(std::array, 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 T prod(matrix::matrix 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 T prod(std::initializer_list init, std::size_t i = 0, std::size_t end = std::string::npos) { std::vector vek; for(auto const& elem : init) vek.push_back(elem); return prod(vek, i, end); } template T prod(std::initializer_list> 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> vek; for(auto const& elem : init) vek.emplace_back(elem); return prod(vek, i, end, j, jend); } template T sum(std::initializer_list init, std::size_t i = 0, std::size_t end = std::string::npos) { std::vector vek; for(auto const& elem : init) vek.push_back(elem); return sum(vek, i, end); } template T sum(std::initializer_list> 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> vek; for(auto const& elem : init) vek.emplace_back(elem); return sum(vek, i, end, j, jend); } } } #endif

### StefanoBelli commented Aug 5, 2019 • edited

I would suggest to use numeric_matrix for MOST usages

```#include "num.h"

using namespace ssynx::num;

int main() {
matrix::numeric_matrix mat { {0} };
matrix::numeric_matrix mat1 { {1,2}, {3,4} };

matrix::numeric_matrix::matrix_type d { {1}, {2} }; //std::vector<std::vector<double>>
matrix::numeric_matrix matd(d);

std::cout << mat + mat1 << mat - mat1 << mat * mat1 << '\n';

auto reduced = matrix::gauss_eliminate(mat);

int row = 0;
int col = 1;
mat.at(row, col);

auto col = mat.col_at(1);
auto row = mat.row_at(0);

auto stdvec = mat.to_vector();

return 0;
}```

to join this conversation on GitHub. Already have an account? Sign in to comment