mat.hpp

浏览该文件的文档.

 
#pragma once
#ifndef XCMATH_MAT_H
#define XCMATH_MAT_H
 
#include <cstddef>
 
#include "../utils/concepts.h"
#include "./declaration.hpp"
#include "./vec.hpp"
 
namespace xcmath {
template <typename _Tp, size_t _rows, size_t _cols = _rows>
class mat : public vec<vec<_Tp, _cols>, _rows> {
   public:
    static constexpr vec2<size_t> shape{_rows, _cols};
 
    using vec<vec<_Tp, _cols>, _rows>::vec;
 
    using vec<vec<_Tp, _cols>, _rows>::operator[];
 
    using vec<vec<_Tp, _cols>, _rows>::begin;
 
    using vec<vec<_Tp, _cols>, _rows>::end;
 
    using vec<vec<_Tp, _cols>, _rows>::operator=;
 
    template <class _T>
    using Self = mat<_T, _rows, _cols>;
 
    mat<_Tp, _rows, _cols>(const vec<vec<_Tp, _cols>, _rows>& o)
        : vec<vec<_Tp, _cols>, _rows>(o) {}
 
    template <size_t _ncols>
    constexpr mat<_Tp, _rows, _ncols> operator^(
        const mat<_Tp, _cols, _ncols>& other) const {
        mat<_Tp, _rows, _ncols> res;
        for (size_t i = 0; i < _rows; i++) {
            for (size_t j = 0; j < _ncols; j++) {
                for (size_t k = 0; k < _cols; k++) {
                    res[i][j] += this->data[i][k] * other[k][j];
                }
            }
        }
        return res;
    }
 
    template <Vec _VTp>
        requires(_VTp::length == _rows) &&
                concepts::Add<typename _VTp::DataType> &&
                concepts::Muitiply<typename _VTp::DataType, _Tp>
    constexpr _VTp operator^(const _VTp& other) const {
        _VTp res;
        for (size_t i = 0; i < _rows; i++) {
            for (size_t j = 0; j < _cols; j++) {
                res[i] += this->data[i][j] * other[j];
            }
        }
        return res;
    }
 
    constexpr mat<_Tp, _cols, _rows> T() const {
        mat<_Tp, _cols, _rows> res;
        for (size_t i = 0; i < _rows; i++) {
            for (size_t j = 0; j < _cols; j++) {
                res[j][i] = this->data[i][j];
            }
        }
        return res;
    }
 
    constexpr mat<_Tp, _cols, _rows> inv() const
        requires(_rows == _cols)
    {
        const _Tp det = this->det();
        mat<_Tp, _rows, _cols> adjugate;
        for (size_t i = 0; i < _rows; ++i) {
            for (size_t j = 0; j < _cols; ++j) {
                // Compute the submatrix (minor)
                mat<_Tp, _rows - 1, _cols - 1> submat;
                size_t sub_i = 0;
                for (size_t k = 0; k < _rows; ++k) {
                    if (k == i) continue;
                    size_t sub_j = 0;
                    for (size_t l = 0; l < _cols; ++l) {
                        if (l == j) continue;
                        submat[sub_i][sub_j++] = this->data[k][l];
                    }
                    ++sub_i;
                }
                // Compute the cofactor (signed minor)
                _Tp cofactor = ((i + j) % 2 == 0 ? 1 : -1) * submat.det();
                adjugate[j][i] = cofactor;  // Transpose for adjugate
            }
        }
        return adjugate / det;
    }
 
    constexpr static mat<_Tp, _rows, _cols> ones() {
        return Self<_Tp>{vec<_Tp, _cols>{_Tp{1}}};
    }
 
    constexpr static mat<_Tp, _rows, _cols> eye() {
        mat<_Tp, _rows, _cols> res;
        size_t n = _rows < _cols ? _rows : _cols;
        for (size_t i = 0; i < n; i++) {
            res[i][i] = 1;
        }
        return res;
    }
 
    constexpr _Tp det() const
        requires(_rows == _cols)
    {
        if constexpr (_rows == 1) {
            return this->data[0][0];
        } else if constexpr (_rows == 2) {
            return this->data[0][0] * this->data[1][1] -
                   this->data[0][1] * this->data[1][0];
        } else {
            _Tp res = 0;
            for (size_t j = 0; j < _rows; j++) {
                mat<_Tp, _rows - 1, _cols - 1> submat;
                for (size_t i = 1; i < _rows; i++) {
                    for (size_t k = 0; k < j; k++) {
                        submat[i - 1][k] = this->data[i][k];
                    }
                    for (size_t k = j + 1; k < _cols; k++) {
                        submat[i - 1][k - 1] = this->data[i][k];
                    }
                }
                res +=
                    ((j % 2 == 0) ? 1 : -1) * this->data[0][j] * submat.det();
            }
            return res;
        }
    }
};
 
}  // namespace xcmath
 
#endif  // XCMATH_MAT_H