matrix-class.hpp

//
// libsemigroups - C++ library for semigroups and monoids
// Copyright (C) 2020-2026 James D. Mitchell
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
//
 
// TODO(1) there're no complete set of init methods for matrices
 
#ifndef LIBSEMIGROUPS_MATRIX_CLASS_HPP_
#define LIBSEMIGROUPS_MATRIX_CLASS_HPP_
 
#include <algorithm>         // for min
#include <array>             // for array
#include <cstddef>           // for size_t
#include <initializer_list>  // for initializer_list
#include <iosfwd>            // for ostringstream
#include <iterator>          // for distance
#include <string>            // for string
#include <tuple>             // for tie
#include <type_traits>       // for false_type, is_signed, true_type
#include <unordered_map>     // for unordered_map
#include <unordered_set>     // for unordered_set
#include <utility>           // for forward, make_pair, pair
#include <vector>            // for vector
 
#include "adapters.hpp"     // for Degree
#include "config.hpp"       // for LIBSEMIGROUPS_PARSED_BY_DOXYGEN
#include "constants.hpp"    // for POSITIVE_INFINITY
#include "debug.hpp"        // for LIBSEMIGROUPS_ASSERT
#include "exception.hpp"    // for LIBSEMIGROUPS_EXCEPTION
#include "is-matrix.hpp"    // for IsMatrix
#include "matrix-view.hpp"  // for RowView
 
#include "detail/containers.hpp"         // for StaticVector1
#include "detail/matrix-common.hpp"      // for entry_repr
#include "detail/matrix-exceptions.hpp"  // for throw_if...
#include "detail/string.hpp"             // for detail::to_string
 
namespace libsemigroups {


  // StaticMatrix with compile time semiring arithmetic

  template <typename PlusOp,
            typename ProdOp,
            typename ZeroOp,
            typename OneOp,
            size_t R,
            size_t C,
            typename Scalar>
  class StaticMatrix
      : public detail::
            MatrixStaticArithmetic<PlusOp, ProdOp, ZeroOp, OneOp, Scalar>,
        public detail::MatrixCommon<
            std::array<Scalar, R * C>,
            StaticMatrix<PlusOp, ProdOp, ZeroOp, OneOp, R, C, Scalar>,
            StaticRowView<PlusOp, ProdOp, ZeroOp, OneOp, C, Scalar>> {
    // StaticMatrix - Aliases - private
 
    using MatrixCommon = ::libsemigroups::detail::MatrixCommon<
        std::array<Scalar, R * C>,
        StaticMatrix<PlusOp, ProdOp, ZeroOp, OneOp, R, C, Scalar>,
        StaticRowView<PlusOp, ProdOp, ZeroOp, OneOp, C, Scalar>>;
    friend MatrixCommon;
 
   public:
    // StaticMatrix - Aliases - public

    using scalar_type = typename MatrixCommon::scalar_type;

    using scalar_reference = typename MatrixCommon::scalar_reference;

    // clang-format off
    // NOLINTNEXTLINE(whitespace/line_length)
    using scalar_const_reference = typename MatrixCommon::scalar_const_reference;
    // clang-format on

    using Row = StaticMatrix<PlusOp, ProdOp, ZeroOp, OneOp, 1, C, Scalar>;

    using RowView = StaticRowView<PlusOp, ProdOp, ZeroOp, OneOp, C, Scalar>;

    using Plus = PlusOp;

    using Prod = ProdOp;

    using Zero = ZeroOp;

    using One = OneOp;

    using iterator = typename MatrixCommon::iterator;

    using const_iterator = typename MatrixCommon::const_iterator;

    static constexpr size_t nr_rows = R;

    static constexpr size_t nr_cols = C;

    // StaticMatrix - Constructors + destructor - public

    explicit StaticMatrix(std::initializer_list<scalar_type> const& c);

    explicit StaticMatrix(
        std::initializer_list<std::initializer_list<scalar_type>> const& m)
        : MatrixCommon(m) {}

    explicit StaticMatrix(std::vector<std::vector<scalar_type>> const& m)
        : MatrixCommon(m) {}

    explicit StaticMatrix(RowView const& rv) : MatrixCommon(rv) {
      static_assert(
          R == 1,
          "cannot construct Matrix with more than one row from RowView!");
    }

    StaticMatrix() = default;

    StaticMatrix(StaticMatrix const&) = default;

    StaticMatrix(StaticMatrix&&) = default;

    StaticMatrix& operator=(StaticMatrix const&) = default;

    StaticMatrix& operator=(StaticMatrix&&) = default;
 
    ~StaticMatrix() = default;
 
#ifndef LIBSEMIGROUPS_PARSED_BY_DOXYGEN
    // For uniformity of interface, the args do nothing
    StaticMatrix(size_t r, size_t c);
 
    // For uniformity of interface, the first arg does nothing
    StaticMatrix(void const* ptr, std::initializer_list<scalar_type> const& c)
        : StaticMatrix(c) {
      (void) ptr;
      LIBSEMIGROUPS_ASSERT(ptr == nullptr);
    }
 
    // For uniformity of interface, the first arg does nothing
    StaticMatrix(
        void const*                                                      ptr,
        std::initializer_list<std::initializer_list<scalar_type>> const& m)
        : StaticMatrix(m) {
      (void) ptr;
      LIBSEMIGROUPS_ASSERT(ptr == nullptr);
    }
 
    // For uniformity of interface, the first arg does nothing
    explicit StaticMatrix(void const* ptr, RowView const& rv)
        : StaticMatrix(rv) {
      (void) ptr;
      LIBSEMIGROUPS_ASSERT(ptr == nullptr);
    }
 
    // For uniformity of interface, no arg used for anything
    StaticMatrix(void const* ptr, size_t r, size_t c) : StaticMatrix(r, c) {
      (void) ptr;
      LIBSEMIGROUPS_ASSERT(ptr == nullptr);
    }
    static StaticMatrix one(size_t n);
    static StaticMatrix one(void const* ptr, size_t n = 0);
#endif  // LIBSEMIGROUPS_PARSED_BY_DOXYGEN

    // StaticMatrix - member function aliases - public
#ifdef LIBSEMIGROUPS_PARSED_BY_DOXYGEN
    scalar_reference operator()(size_t r, size_t c);

    scalar_reference at(size_t r, size_t c);

    scalar_const_reference operator()(size_t r, size_t c) const;

    scalar_const_reference at(size_t r, size_t c) const;

    iterator begin() noexcept;

    iterator end();

    const_iterator cbegin() const noexcept;

    const_iterator cend();

    bool operator==(StaticMatrix const& that) const;

    bool operator==(RowView const& that) const;

    bool operator!=(StaticMatrix const& that) const;

    bool operator!=(RowView const& that) const;

    bool operator<(StaticMatrix const& that) const;

    bool operator<(RowView const& that) const;

    bool operator>(StaticMatrix const& that) const;

    std::pair<scalar_type, scalar_type> coords(const_iterator it) const;

    size_t number_of_rows() const noexcept;

    size_t number_of_cols() const noexcept;

    StaticMatrix plus_no_checks(StaticMatrix const& y) const;

    StaticMatrix operator+(StaticMatrix const& that);

    void plus_inplace_no_checks(StaticMatrix const& that);

    void plus_inplace_no_checks(RowView const& that);

    void operator+=(StaticMatrix const& that);

    void operator+=(RowView const& that);

    void operator+=(scalar_type a);

    StaticMatrix product_no_checks(StaticMatrix const& that);

    StaticMatrix operator*(StaticMatrix const& that);

    void operator*=(scalar_type a);

    void product_inplace_no_checks(StaticMatrix const& x,
                                   StaticMatrix const& y);

    void product_inplace(StaticMatrix const& x, StaticMatrix const& y);

    RowView row_no_checks(size_t i) const;

    RowView row(size_t i) const;

    template <typename T>
    void rows(T& x) const;

    void swap(StaticMatrix& that) noexcept;

    void transpose();

    void transpose_no_checks();

    static StaticMatrix one() const;

    size_t hash_value() const;

    template <typename U>
    bool operator<=(U const& that) const;

    template <typename U>
    bool operator>=(U const& that) const;

    StaticMatrix operator*(scalar_type a);

    StaticMatrix operator+(scalar_type a);

    scalar_type scalar_zero() const noexcept;

    scalar_type scalar_one() const noexcept;

    semiring_type const* semiring() const noexcept;
 
#else
    using MatrixCommon::at;
    using MatrixCommon::begin;
    using MatrixCommon::cbegin;
    using MatrixCommon::cend;
    using MatrixCommon::coords;
    using MatrixCommon::end;
    using MatrixCommon::hash_value;
    using MatrixCommon::number_of_cols;
    using MatrixCommon::number_of_rows;
    using MatrixCommon::one;
    using MatrixCommon::operator!=;
    using MatrixCommon::operator();
    using MatrixCommon::operator*;
    using MatrixCommon::operator*=;
    using MatrixCommon::operator+;
    using MatrixCommon::operator+=;
    using MatrixCommon::operator<;  // NOLINT(whitespace/operators)
    using MatrixCommon::operator<=;
    using MatrixCommon::operator==;
    using MatrixCommon::operator>;  // NOLINT(whitespace/operators)
    using MatrixCommon::operator>=;
    using MatrixCommon::product_inplace_no_checks;
    using MatrixCommon::row;
    using MatrixCommon::row_no_checks;
    using MatrixCommon::rows;
    using MatrixCommon::scalar_one;
    using MatrixCommon::scalar_zero;
    using MatrixCommon::semiring;
    using MatrixCommon::swap;
    using MatrixCommon::transpose;
    using MatrixCommon::transpose_no_checks;
#endif  // LIBSEMIGROUPS_PARSED_BY_DOXYGEN
 
   private:
    // StaticMatrix - implementation of MatrixCommon requirements - private
 
    static constexpr size_t number_of_rows_impl() noexcept {
      return R;
    }
    static constexpr size_t number_of_cols_impl() noexcept {
      return C;
    }
  };  // class StaticMatrix

  // DynamicMatrix with compile time semiring arithmetic

  template <typename PlusOp,
            typename ProdOp,
            typename ZeroOp,
            typename OneOp,
            typename Scalar>
  class DynamicMatrix<PlusOp, ProdOp, ZeroOp, OneOp, Scalar>
      : public detail::MatrixDynamicDim<Scalar>,
        public detail::MatrixCommon<
            std::vector<Scalar>,
            DynamicMatrix<PlusOp, ProdOp, ZeroOp, OneOp, Scalar>,
            DynamicRowView<PlusOp, ProdOp, ZeroOp, OneOp, Scalar>>,
        public detail::
            MatrixStaticArithmetic<PlusOp, ProdOp, ZeroOp, OneOp, Scalar> {
    using MatrixDynamicDim = ::libsemigroups::detail::MatrixDynamicDim<Scalar>;
    using MatrixCommon     = ::libsemigroups::detail::MatrixCommon<
        std::vector<Scalar>,
        DynamicMatrix<PlusOp, ProdOp, ZeroOp, OneOp, Scalar>,
        DynamicRowView<PlusOp, ProdOp, ZeroOp, OneOp, Scalar>>;
    friend MatrixCommon;
 
   public:
    using scalar_type = typename MatrixCommon::scalar_type;

    using scalar_reference = typename MatrixCommon::scalar_reference;

    // clang-format off
    // NOLINTNEXTLINE(whitespace/line_length)
    using scalar_const_reference = typename MatrixCommon::scalar_const_reference;
    // clang-format on

    using Row = DynamicMatrix;

    using RowView = DynamicRowView<PlusOp, ProdOp, ZeroOp, OneOp, Scalar>;

    using Plus = PlusOp;

    using Prod = ProdOp;

    using Zero = ZeroOp;

    using One = OneOp;

    using semiring_type = void;

    DynamicMatrix() = default;

    DynamicMatrix(DynamicMatrix const&) = default;

    DynamicMatrix(DynamicMatrix&&) = default;

    DynamicMatrix& operator=(DynamicMatrix const&) = default;

    DynamicMatrix& operator=(DynamicMatrix&&) = default;
 
    ~DynamicMatrix() = default;

    DynamicMatrix(size_t r, size_t c) : MatrixDynamicDim(r, c), MatrixCommon() {
      resize(number_of_rows(), number_of_cols());
    }

    explicit DynamicMatrix(std::initializer_list<scalar_type> const& c)
        : MatrixDynamicDim(1, c.size()), MatrixCommon(c) {}

    explicit DynamicMatrix(
        std::initializer_list<std::initializer_list<scalar_type>> const& m)
        : MatrixDynamicDim(m.size(), std::empty(m) ? 0 : m.begin()->size()),
          MatrixCommon(m) {}

    explicit DynamicMatrix(std::vector<std::vector<scalar_type>> const& m)
        : MatrixDynamicDim(m.size(), std::empty(m) ? 0 : m.begin()->size()),
          MatrixCommon(m) {}

    explicit DynamicMatrix(RowView const& rv)
        : MatrixDynamicDim(1, rv.size()), MatrixCommon(rv) {}
 
#ifndef LIBSEMIGROUPS_PARSED_BY_DOXYGEN
    // For uniformity of interface, the first arg does nothing
    DynamicMatrix(void const* ptr, size_t r, size_t c) : DynamicMatrix(r, c) {
      (void) ptr;
      LIBSEMIGROUPS_ASSERT(ptr == nullptr);
    }
 
    // For uniformity of interface, the first arg does nothing
    DynamicMatrix(void const* ptr, std::initializer_list<scalar_type> const& c)
        : DynamicMatrix(c) {
      (void) ptr;
      LIBSEMIGROUPS_ASSERT(ptr == nullptr);
    }
 
    // For uniformity of interface, the first arg does nothing
    DynamicMatrix(
        void const*                                                      ptr,
        std::initializer_list<std::initializer_list<scalar_type>> const& m)
        : DynamicMatrix(m) {
      (void) ptr;
      LIBSEMIGROUPS_ASSERT(ptr == nullptr);
    }
 
    static DynamicMatrix one(void const* ptr, size_t n) {
      (void) ptr;
      LIBSEMIGROUPS_ASSERT(ptr == nullptr);
      return one(n);
    }
#endif  // LIBSEMIGROUPS_PARSED_BY_DOXYGEN

    static DynamicMatrix one(size_t n) {
      DynamicMatrix x(n, n);
      std::fill(x.begin(), x.end(), ZeroOp()());
      for (size_t r = 0; r < n; ++r) {
        x(r, r) = OneOp()();
      }
      return x;
    }
 
#ifdef LIBSEMIGROUPS_PARSED_BY_DOXYGEN
    scalar_reference at(size_t r, size_t c);

    scalar_reference at(size_t r, size_t c) const;

    iterator begin() noexcept;

    const_iterator cbegin() noexcept;

    const_iterator cend() noexcept;

    std::pair<scalar_type, scalar_type> coords(const_iterator it) const;

    iterator end() noexcept;

    size_t hash_value() const;

    size_t number_of_cols() const noexcept;

    size_t number_of_rows() const noexcept;

    bool operator!=(DynamicMatrix const& that) const;

    bool operator!=(RowView const& that) const;

    scalar_reference operator()(size_t r, size_t c);

    scalar_const_reference operator()(size_t r, size_t c) const;
    DynamicMatrix operator*(scalar_type a);

    DynamicMatrix product_no_checks(DynamicMatrix const& that);

    DynamicMatrix operator*(DynamicMatrix const& that);

    void operator*=(scalar_type a);

    DynamicMatrix plus_no_checks(DynamicMatrix const& y) const;

    DynamicMatrix operator+(DynamicMatrix const& that);

    void plus_inplace_no_checks(DynamicMatrix const& that);

    void plus_inplace_no_checks(RowView const& that);

    void operator+=(DynamicMatrix const& that);

    void operator+=(RowView const& that);

    void operator+=(scalar_type a);

    bool operator<(DynamicMatrix const& that) const;

    bool operator<(RowView const& that) const;

    template <typename T>
    bool operator<=(T const& that) const;

    bool operator==(DynamicMatrix const& that) const;

    bool operator==(RowView const& that) const;

    bool operator>(DynamicMatrix const& that) const;

    template <typename T>
    bool operator>=(T const& that) const;

    void product_inplace_no_checks(DynamicMatrix const& x,
                                   DynamicMatrix const& y);

    void product_inplace(DynamicMatrix const& x, DynamicMatrix const& y);

    RowView row(size_t i) const;

    RowView row_no_checks(size_t i) const;

    template <typename T>
    void rows(T& x) const;

    scalar_type scalar_one() const noexcept;

    scalar_type scalar_zero() const noexcept;

    semiring_type const* semiring() const noexcept;

    void transpose();

    void transpose_no_checks();
#else
    using MatrixCommon::at;
    using MatrixCommon::begin;
    using MatrixCommon::cbegin;
    using MatrixCommon::cend;
    using MatrixCommon::coords;
    using MatrixCommon::end;
    using MatrixCommon::hash_value;
    using MatrixCommon::number_of_cols;
    using MatrixCommon::number_of_rows;
    using MatrixCommon::one;
    using MatrixCommon::operator!=;
    using MatrixCommon::operator();
    using MatrixCommon::operator*;
    using MatrixCommon::operator*=;
    using MatrixCommon::operator+;
    using MatrixCommon::operator+=;
    using MatrixCommon::operator<;  // NOLINT(whitespace/operators)
    using MatrixCommon::operator<=;
    using MatrixCommon::operator==;
    using MatrixCommon::operator>;  // NOLINT(whitespace/operators)
    using MatrixCommon::operator>=;
    using MatrixCommon::product_inplace_no_checks;
    using MatrixCommon::row;
    using MatrixCommon::row_no_checks;
    using MatrixCommon::rows;
    using MatrixCommon::scalar_one;
    using MatrixCommon::scalar_zero;
    using MatrixCommon::semiring;
    // using MatrixCommon::swap; // Don't want this use the one below.
    using MatrixCommon::transpose;
    using MatrixCommon::transpose_no_checks;
#endif  // LIBSEMIGROUPS_PARSED_BY_DOXYGEN

    void swap(DynamicMatrix& that) noexcept {
      static_cast<MatrixDynamicDim&>(*this).swap(
          static_cast<MatrixDynamicDim&>(that));
      static_cast<MatrixCommon&>(*this).swap(static_cast<MatrixCommon&>(that));
    }
 
   private:
    using MatrixCommon::resize;
  };

  // DynamicMatrix with runtime semiring arithmetic

  template <typename Semiring, typename Scalar>
  class DynamicMatrix<Semiring, Scalar>
      : public detail::MatrixDynamicDim<Scalar>,
        public detail::MatrixCommon<std::vector<Scalar>,
                                    DynamicMatrix<Semiring, Scalar>,
                                    DynamicRowView<Semiring, Scalar>,
                                    Semiring> {
    using MatrixCommon = detail::MatrixCommon<std::vector<Scalar>,
                                              DynamicMatrix<Semiring, Scalar>,
                                              DynamicRowView<Semiring, Scalar>,
                                              Semiring>;
    friend MatrixCommon;
    using MatrixDynamicDim = ::libsemigroups::detail::MatrixDynamicDim<Scalar>;
 
   public:
    using scalar_type = typename MatrixCommon::scalar_type;

    using scalar_reference = typename MatrixCommon::scalar_reference;

    // clang-format off
    // NOLINTNEXTLINE(whitespace/line_length)
    using scalar_const_reference = typename MatrixCommon::scalar_const_reference;
    // clang-format on

    using Row = DynamicMatrix;

    using RowView = DynamicRowView<Semiring, Scalar>;
 
    friend RowView;

    using semiring_type = Semiring;

    DynamicMatrix() = delete;

    DynamicMatrix(DynamicMatrix const&) = default;

    DynamicMatrix(DynamicMatrix&&) = default;

    DynamicMatrix& operator=(DynamicMatrix const&) = default;

    DynamicMatrix& operator=(DynamicMatrix&&) = default;

    DynamicMatrix(Semiring const* sr, size_t r, size_t c)
        : MatrixDynamicDim(r, c), MatrixCommon(), _semiring(sr) {
      resize(number_of_rows(), number_of_cols());
    }

    DynamicMatrix(
        Semiring const*                                                  sr,
        std::initializer_list<std::initializer_list<scalar_type>> const& rws);

    explicit DynamicMatrix(Semiring const*                              sr,
                           std::vector<std::vector<scalar_type>> const& rws)
        : MatrixDynamicDim(rws.size(), (rws.empty() ? 0 : rws.begin()->size())),
          MatrixCommon(rws),
          _semiring(sr) {}

    explicit DynamicMatrix(Semiring const*                           sr,
                           std::initializer_list<scalar_type> const& rw)
        : MatrixDynamicDim(1, rw.size()), MatrixCommon(rw), _semiring(sr) {}

    explicit DynamicMatrix(RowView const& rv)
        : MatrixDynamicDim(1, rv.size()),
          MatrixCommon(rv),
          _semiring(rv._matrix->semiring()) {}

    // No static DynamicMatrix::one(size_t n) because we need a semiring!
    static DynamicMatrix one(Semiring const* semiring, size_t n);
 
    ~DynamicMatrix() = default;
 
#ifdef LIBSEMIGROUPS_PARSED_BY_DOXYGEN
    scalar_reference at(size_t r, size_t c);

    scalar_reference at(size_t r, size_t c) const;

    iterator begin() noexcept;

    const_iterator cbegin() noexcept;

    const_iterator cend() noexcept;

    std::pair<scalar_type, scalar_type> coords(const_iterator it) const;

    iterator end() noexcept;

    size_t hash_value() const;

    size_t number_of_cols() const noexcept;

    size_t number_of_rows() const noexcept;

    bool operator!=(DynamicMatrix const& that) const;

    bool operator!=(RowView const& that) const;

    scalar_reference operator()(size_t r, size_t c);

    scalar_const_reference operator()(size_t r, size_t c) const;

    DynamicMatrix operator*(scalar_type a);

    DynamicMatrix product_no_checks(DynamicMatrix const& that);

    DynamicMatrix operator*(DynamicMatrix const& that);

    void operator*=(scalar_type a);

    DynamicMatrix operator+(DynamicMatrix const& that);

    void plus_inplace_no_checks(DynamicMatrix const& that);

    void plus_inplace_no_checks(RowView const& that);

    void operator+=(DynamicMatrix const& that);

    void operator+=(RowView const& that);

    void operator+=(scalar_type a);

    bool operator<(DynamicMatrix const& that) const;

    bool operator<(RowView const& that) const;

    template <typename T>
    bool operator<=(T const& that) const;

    bool operator==(DynamicMatrix const& that) const;

    bool operator==(RowView const& that) const;

    bool operator>(DynamicMatrix const& that) const;

    template <typename T>
    bool operator>=(T const& that) const;

    void product_inplace_no_checks(DynamicMatrix const& x,
                                   DynamicMatrix const& y);

    void product_inplace(DynamicMatrix const& x, DynamicMatrix const& y);

    RowView row(size_t i) const;

    RowView row_no_checks(size_t i) const;

    template <typename T>
    void rows(T& x) const;

    scalar_type scalar_one() const noexcept;

    scalar_type scalar_zero() const noexcept;

    semiring_type const* semiring() const noexcept;

    void transpose();

    void transpose_no_checks();
#else
    using MatrixCommon::at;
    using MatrixCommon::begin;
    using MatrixCommon::cbegin;
    using MatrixCommon::cend;
    using MatrixCommon::coords;
    using MatrixCommon::end;
    using MatrixCommon::hash_value;
    using MatrixCommon::number_of_cols;
    using MatrixCommon::number_of_rows;
    using MatrixCommon::one;
    using MatrixCommon::operator!=;
    using MatrixCommon::operator();
    using MatrixCommon::operator*;
    using MatrixCommon::operator*=;
    using MatrixCommon::operator+;
    using MatrixCommon::operator+=;
    using MatrixCommon::operator<;  // NOLINT(whitespace/operators)
    using MatrixCommon::operator<=;
    using MatrixCommon::operator==;
    using MatrixCommon::operator>;  // NOLINT(whitespace/operators)
    using MatrixCommon::operator>=;
    using MatrixCommon::product_inplace_no_checks;
    using MatrixCommon::row;
    using MatrixCommon::row_no_checks;
    using MatrixCommon::rows;
    using MatrixCommon::scalar_one;
    using MatrixCommon::scalar_zero;
    using MatrixCommon::semiring;
    // using MatrixCommon::swap; // Don't want this use the one below.
    using MatrixCommon::transpose;
    using MatrixCommon::transpose_no_checks;
#endif  // LIBSEMIGROUPS_PARSED_BY_DOXYGEN

    void swap(DynamicMatrix& that) noexcept;
 
   private:
    using MatrixCommon::resize;
 
    scalar_type plus_no_checks_impl(scalar_type x,
                                    scalar_type y) const noexcept {
      return _semiring->plus_no_checks(x, y);
    }
 
    scalar_type product_no_checks_impl(scalar_type x,
                                       scalar_type y) const noexcept {
      return _semiring->product_no_checks(x, y);
    }
 
    scalar_type one_impl() const noexcept {
      return _semiring->scalar_one();
    }
 
    scalar_type zero_impl() const noexcept {
      return _semiring->scalar_zero();
    }
 
    Semiring const* semiring_impl() const noexcept {
      return _semiring;
    }
 
    Semiring const* _semiring;
  };

  // Helper structs to check if matrix is static, or has a pointer to a
  // semiring
 
  namespace detail {
    template <typename PlusOp,
              typename ProdOp,
              typename ZeroOp,
              typename OneOp,
              size_t R,
              size_t C,
              typename Scalar>
    struct IsStaticMatrixHelper<
        StaticMatrix<PlusOp, ProdOp, ZeroOp, OneOp, R, C, Scalar>>
        : std::true_type {};
 
    template <typename Semiring, typename Scalar>
    struct IsMatWithSemiringHelper<DynamicMatrix<Semiring, Scalar>>
        : std::true_type {};
  }  // namespace detail

  // Boolean matrices - compile time semiring arithmetic


  struct BooleanPlus {
    constexpr bool operator()(bool x, bool y) const noexcept {
      return x || y;
    }
  };

  struct BooleanProd {
    constexpr bool operator()(bool x, bool y) const noexcept {
      return x & y;
    }
  };

  struct BooleanOne {
    constexpr bool operator()() const noexcept {
      return true;
    }
  };

  struct BooleanZero {
    constexpr bool operator()() const noexcept {
      return false;
    }
  };

  // The use of `int` rather than `bool` as the scalar type for dynamic
  // boolean matrices is intentional, because the bit iterators implemented in
  // std::vector<bool> entail a significant performance penalty.
  using DynamicBMat
      = DynamicMatrix<BooleanPlus, BooleanProd, BooleanZero, BooleanOne, int>;

  template <size_t R, size_t C>
  using StaticBMat = StaticMatrix<BooleanPlus,
                                  BooleanProd,
                                  BooleanZero,
                                  BooleanOne,
                                  R,
                                  C,
                                  int>;

  // FLS + JDM considered adding BMat8 and decided it wasn't a good idea.
  template <size_t R = 0, size_t C = R>
  using BMat
      = std::conditional_t<R == 0 || C == 0, DynamicBMat, StaticBMat<R, C>>;
 
  namespace detail {
    template <size_t R, size_t C>
    struct IsBMatHelper<StaticBMat<R, C>> : std::true_type {};
 
    template <>
    struct IsBMatHelper<DynamicBMat> : std::true_type {};
  }  // namespace detail

  // Integer matrices - compile time semiring arithmetic


  //! \tparam Scalar the type of the entries in the matrix.
  template <typename Scalar>
  struct IntegerPlus {
    //! \exceptions
    //! \noexcept
    constexpr Scalar operator()(Scalar x, Scalar y) const noexcept {
      return x + y;
    }
  };

  //! \tparam Scalar the type of the entries in the matrix.
  template <typename Scalar>
  struct IntegerProd {
    //! \exceptions
    //! \noexcept
    constexpr Scalar operator()(Scalar x, Scalar y) const noexcept {
      return x * y;
    }
  };

  //! the additive identity of the integer semiring.
  template <typename Scalar>
  struct IntegerZero {
    //! \exceptions
    //! \noexcept
    constexpr Scalar operator()() const noexcept {
      return 0;
    }
  };

  //! the multiplicative identity of the integer semiring.
  template <typename Scalar>
  struct IntegerOne {
    //! \exceptions
    //! \noexcept
    constexpr Scalar operator()() const noexcept {
      return 1;
    }
  };

  template <typename Scalar>
  using DynamicIntMat = DynamicMatrix<IntegerPlus<Scalar>,
                                      IntegerProd<Scalar>,
                                      IntegerZero<Scalar>,
                                      IntegerOne<Scalar>,
                                      Scalar>;

  template <size_t R, size_t C, typename Scalar>
  using StaticIntMat = StaticMatrix<IntegerPlus<Scalar>,
                                    IntegerProd<Scalar>,
                                    IntegerZero<Scalar>,
                                    IntegerOne<Scalar>,
                                    R,
                                    C,
                                    Scalar>;

  template <size_t R = 0, size_t C = R, typename Scalar = int>
  using IntMat = std::conditional_t<R == 0 || C == 0,
                                    DynamicIntMat<Scalar>,
                                    StaticIntMat<R, C, Scalar>>;
  namespace detail {
    template <size_t R, size_t C, typename Scalar>
    struct IsIntMatHelper<StaticIntMat<R, C, Scalar>> : std::true_type {};
 
    template <typename Scalar>
    struct IsIntMatHelper<DynamicIntMat<Scalar>> : std::true_type {};
  }  // namespace detail

  // Max-plus matrices

  // Static arithmetic
  template <typename Scalar>
  struct MaxPlusPlus {
    static_assert(std::is_signed<Scalar>::value,
                  "MaxPlus requires a signed integer type as parameter!");
    Scalar operator()(Scalar x, Scalar y) const noexcept;
  };

  //! integer type).
  template <typename Scalar>
  struct MaxPlusProd {
    static_assert(std::is_signed<Scalar>::value,
                  "MaxPlus requires a signed integer type as parameter!");
    Scalar operator()(Scalar x, Scalar y) const noexcept;
  };

  //! integer type).
  template <typename Scalar>
  struct MaxPlusZero {
    static_assert(std::is_signed<Scalar>::value,
                  "MaxPlus requires a signed integer type as parameter!");
    //! \exceptions
    //! \noexcept
    constexpr Scalar operator()() const noexcept {
      return NEGATIVE_INFINITY;
    }
  };

  template <typename Scalar>
  using DynamicMaxPlusMat = DynamicMatrix<MaxPlusPlus<Scalar>,
                                          MaxPlusProd<Scalar>,
                                          MaxPlusZero<Scalar>,
                                          IntegerZero<Scalar>,
                                          Scalar>;

  template <size_t R, size_t C, typename Scalar>
  using StaticMaxPlusMat = StaticMatrix<MaxPlusPlus<Scalar>,
                                        MaxPlusProd<Scalar>,
                                        MaxPlusZero<Scalar>,
                                        IntegerZero<Scalar>,
                                        R,
                                        C,
                                        Scalar>;

  template <size_t R = 0, size_t C = R, typename Scalar = int>
  using MaxPlusMat = std::conditional_t<R == 0 || C == 0,
                                        DynamicMaxPlusMat<Scalar>,
                                        StaticMaxPlusMat<R, C, Scalar>>;
 
  namespace detail {
    template <size_t R, size_t C, typename Scalar>
    struct IsMaxPlusMatHelper<StaticMaxPlusMat<R, C, Scalar>> : std::true_type {
    };
 
    template <typename Scalar>
    struct IsMaxPlusMatHelper<DynamicMaxPlusMat<Scalar>> : std::true_type {};
  }  // namespace detail

  // Min-plus matrices


  //!
  template <typename Scalar>
  struct MinPlusPlus {
    static_assert(std::is_signed<Scalar>::value,
                  "MinPlus requires a signed integer type as parameter!");
    Scalar operator()(Scalar x, Scalar y) const noexcept;
  };

  //! integer type).
  template <typename Scalar>
  struct MinPlusProd {
    static_assert(std::is_signed<Scalar>::value,
                  "MinPlus requires a signed integer type as parameter!");
    Scalar operator()(Scalar x, Scalar y) const noexcept;
  };

  //! integer type).
  template <typename Scalar>
  struct MinPlusZero {
    static_assert(std::is_signed<Scalar>::value,
                  "MinPlus requires a signed integer type as parameter!");
    //! \exceptions
    //! \noexcept
    constexpr Scalar operator()() const noexcept {
      return POSITIVE_INFINITY;
    }
  };

  template <typename Scalar>
  using DynamicMinPlusMat = DynamicMatrix<MinPlusPlus<Scalar>,
                                          MinPlusProd<Scalar>,
                                          MinPlusZero<Scalar>,
                                          IntegerZero<Scalar>,
                                          Scalar>;

  template <size_t R, size_t C, typename Scalar>
  using StaticMinPlusMat = StaticMatrix<MinPlusPlus<Scalar>,
                                        MinPlusProd<Scalar>,
                                        MinPlusZero<Scalar>,
                                        IntegerZero<Scalar>,
                                        R,
                                        C,
                                        Scalar>;
  template <size_t R = 0, size_t C = R, typename Scalar = int>
  using MinPlusMat = std::conditional_t<R == 0 || C == 0,
                                        DynamicMinPlusMat<Scalar>,
                                        StaticMinPlusMat<R, C, Scalar>>;
 
  namespace detail {
    template <size_t R, size_t C, typename Scalar>
    struct IsMinPlusMatHelper<StaticMinPlusMat<R, C, Scalar>> : std::true_type {
    };
 
    template <typename Scalar>
    struct IsMinPlusMatHelper<DynamicMinPlusMat<Scalar>> : std::true_type {};
  }  // namespace detail

  // Max-plus matrices with threshold


  //! integer type).
  template <size_t T, typename Scalar>
  struct MaxPlusTruncProd {
    static_assert(std::is_signed<Scalar>::value,
                  "MaxPlus requires a signed integer type as parameter!");
    Scalar operator()(Scalar x, Scalar y) const noexcept;
  };

  //! signed integer type (defaults to \c int).
  template <typename Scalar = int>
  class MaxPlusTruncSemiring {
    static_assert(std::is_signed<Scalar>::value,
                  "MaxPlus requires a signed integer type as parameter!");
 
   public:
    MaxPlusTruncSemiring() = delete;

    MaxPlusTruncSemiring(MaxPlusTruncSemiring const&) noexcept = default;

    MaxPlusTruncSemiring(MaxPlusTruncSemiring&&) noexcept = default;

    MaxPlusTruncSemiring& operator=(MaxPlusTruncSemiring const&) noexcept
        = default;

    MaxPlusTruncSemiring& operator=(MaxPlusTruncSemiring&&) noexcept = default;
 
    ~MaxPlusTruncSemiring() = default;

    explicit MaxPlusTruncSemiring(Scalar threshold);

    //! \exceptions
    //! \noexcept
    static constexpr Scalar scalar_one() noexcept {
      return 0;
    }

    //! \exceptions
    //! \noexcept
    static constexpr Scalar scalar_zero() noexcept {
      return NEGATIVE_INFINITY;
    }

    Scalar product_no_checks(Scalar x, Scalar y) const noexcept;

    Scalar plus_no_checks(Scalar x, Scalar y) const noexcept;

    //! \complexity
    //! Constant.
    Scalar threshold() const noexcept {
      return _threshold;
    }
 
   public:
    Scalar const _threshold;
  };

  template <size_t T, typename Scalar>
  using DynamicMaxPlusTruncMat = DynamicMatrix<MaxPlusPlus<Scalar>,
                                               MaxPlusTruncProd<T, Scalar>,
                                               MaxPlusZero<Scalar>,
                                               IntegerZero<Scalar>,
                                               Scalar>;

  template <size_t T, size_t R, size_t C, typename Scalar>
  using StaticMaxPlusTruncMat = StaticMatrix<MaxPlusPlus<Scalar>,
                                             MaxPlusTruncProd<T, Scalar>,
                                             MaxPlusZero<Scalar>,
                                             IntegerZero<Scalar>,
                                             R,
                                             C,
                                             Scalar>;
  template <size_t T = 0, size_t R = 0, size_t C = R, typename Scalar = int>
  using MaxPlusTruncMat = std::conditional_t<
      R == 0 || C == 0,
      std::conditional_t<T == 0,
                         DynamicMatrix<MaxPlusTruncSemiring<Scalar>, Scalar>,
                         DynamicMaxPlusTruncMat<T, Scalar>>,
      StaticMaxPlusTruncMat<T, R, C, Scalar>>;
 
  namespace detail {
    template <size_t T, size_t R, size_t C, typename Scalar>
    struct IsMaxPlusTruncMatHelper<StaticMaxPlusTruncMat<T, R, C, Scalar>>
        : std::true_type {
      static constexpr Scalar threshold = T;
    };
 
    template <size_t T, typename Scalar>
    struct IsMaxPlusTruncMatHelper<DynamicMaxPlusTruncMat<T, Scalar>>
        : std::true_type {
      static constexpr Scalar threshold = T;
    };
 
    template <typename Scalar>
    struct IsMaxPlusTruncMatHelper<
        DynamicMatrix<MaxPlusTruncSemiring<Scalar>, Scalar>> : std::true_type {
      static constexpr Scalar threshold = UNDEFINED;
    };
 
    template <typename T>
    struct IsTruncMatHelper<T, std::enable_if_t<IsMaxPlusTruncMat<T>>>
        : std::true_type {
      static constexpr typename T::scalar_type threshold
          = IsMaxPlusTruncMatHelper<T>::threshold;
    };
  }  // namespace detail

  // Min-plus matrices with threshold


  //! \tparam Scalar the type of the values in the semiring.
  template <size_t T, typename Scalar>
  struct MinPlusTruncProd {
    Scalar operator()(Scalar x, Scalar y) const noexcept;
  };  // MinPlusTruncProd

  //! integral type.
  template <typename Scalar = int>
  class MinPlusTruncSemiring {
    static_assert(std::is_integral<Scalar>::value,
                  "MinPlus requires an integral type as parameter!");
 
   public:
    MinPlusTruncSemiring() = delete;

    MinPlusTruncSemiring(MinPlusTruncSemiring const&) noexcept = default;

    MinPlusTruncSemiring(MinPlusTruncSemiring&&) noexcept = default;

    MinPlusTruncSemiring& operator=(MinPlusTruncSemiring const&) noexcept
        = default;

    MinPlusTruncSemiring& operator=(MinPlusTruncSemiring&&) noexcept = default;

    explicit MinPlusTruncSemiring(Scalar threshold);

    //! \exceptions
    //! \noexcept
    static constexpr Scalar scalar_one() noexcept {
      return 0;
    }

    //! \noexcept
    // TODO(1) These mem fns (one and zero) aren't needed?
    static constexpr Scalar scalar_zero() noexcept {
      return POSITIVE_INFINITY;
    }

    Scalar product_no_checks(Scalar x, Scalar y) const noexcept;

    Scalar plus_no_checks(Scalar x, Scalar y) const noexcept;

    //! \complexity
    //! Constant.
    Scalar threshold() const noexcept {
      return _threshold;
    }
 
   public:
    Scalar const _threshold;
  };

  template <size_t T, typename Scalar>
  using DynamicMinPlusTruncMat = DynamicMatrix<MinPlusPlus<Scalar>,
                                               MinPlusTruncProd<T, Scalar>,
                                               MinPlusZero<Scalar>,
                                               IntegerZero<Scalar>,
                                               Scalar>;

  template <size_t T, size_t R, size_t C, typename Scalar>
  using StaticMinPlusTruncMat = StaticMatrix<MinPlusPlus<Scalar>,
                                             MinPlusTruncProd<T, Scalar>,
                                             MinPlusZero<Scalar>,
                                             IntegerZero<Scalar>,
                                             R,
                                             C,
                                             Scalar>;

  template <size_t T = 0, size_t R = 0, size_t C = R, typename Scalar = int>
  using MinPlusTruncMat = std::conditional_t<
      R == 0 || C == 0,
      std::conditional_t<T == 0,
                         DynamicMatrix<MinPlusTruncSemiring<Scalar>, Scalar>,
                         DynamicMinPlusTruncMat<T, Scalar>>,
      StaticMinPlusTruncMat<T, R, C, Scalar>>;
 
  namespace detail {
    template <size_t T, size_t R, size_t C, typename Scalar>
    struct IsMinPlusTruncMatHelper<StaticMinPlusTruncMat<T, R, C, Scalar>>
        : std::true_type {
      static constexpr Scalar threshold = T;
    };
 
    template <size_t T, typename Scalar>
    struct IsMinPlusTruncMatHelper<DynamicMinPlusTruncMat<T, Scalar>>
        : std::true_type {
      static constexpr Scalar threshold = T;
    };
 
    template <typename Scalar>
    struct IsMinPlusTruncMatHelper<
        DynamicMatrix<MinPlusTruncSemiring<Scalar>, Scalar>> : std::true_type {
      static constexpr Scalar threshold = UNDEFINED;
    };
 
    template <typename T>
    struct IsTruncMatHelper<T, std::enable_if_t<IsMinPlusTruncMat<T>>>
        : std::true_type {
      static constexpr typename T::scalar_type threshold
          = IsMinPlusTruncMatHelper<T>::threshold;
    };
  }  // namespace detail

  // NTP matrices

 
  namespace detail {
    template <size_t T, size_t P, typename Scalar>
    constexpr Scalar thresholdperiod(Scalar x) noexcept;
 
    template <typename Scalar>
    inline Scalar thresholdperiod(Scalar x,
                                  Scalar threshold,
                                  Scalar period) noexcept;
  }  // namespace detail

  // Static arithmetic
  template <size_t T, size_t P, typename Scalar>
  struct NTPPlus {
    //! \exceptions
    //! \noexcept
    constexpr Scalar operator()(Scalar x, Scalar y) const noexcept {
      return detail::thresholdperiod<T, P>(x + y);
    }
  };

  //! \tparam Scalar the type of the values in the semiring.
  template <size_t T, size_t P, typename Scalar>
  struct NTPProd {
    //! \exceptions
    //! \noexcept
    constexpr Scalar operator()(Scalar x, Scalar y) const noexcept {
      return detail::thresholdperiod<T, P>(x * y);
    }
  };

  // Dynamic arithmetic
  template <typename Scalar = size_t>
  class NTPSemiring {
   public:
    // Deleted to avoid uninitialised values of period and threshold.
    NTPSemiring() = delete;

    NTPSemiring(NTPSemiring const&) = default;

    NTPSemiring(NTPSemiring&&) = default;

    NTPSemiring& operator=(NTPSemiring const&) = default;

    NTPSemiring& operator=(NTPSemiring&&) = default;
 
    ~NTPSemiring() = default;

    //! \complexity
    //! Constant.
    NTPSemiring(Scalar t, Scalar p) : _period(p), _threshold(t) {
      if constexpr (std::is_signed<Scalar>::value) {
        if (t < 0) {
          LIBSEMIGROUPS_EXCEPTION(
              "expected non-negative value for 1st argument, found {}", t);
        }
      }
      if (p <= 0) {
        LIBSEMIGROUPS_EXCEPTION(
            "expected positive value for 2nd argument, found {}", p);
      }
    }

    //! \exceptions
    //! \noexcept
    static constexpr Scalar scalar_one() noexcept {
      return 1;
    }

    //! \complexity
    //! Constant.
    static constexpr Scalar scalar_zero() noexcept {
      return 0;
    }

    //! \complexity
    //! Constant.
    Scalar product_no_checks(Scalar x, Scalar y) const noexcept {
      LIBSEMIGROUPS_ASSERT(x >= 0 && x <= _period + _threshold - 1);
      LIBSEMIGROUPS_ASSERT(y >= 0 && y <= _period + _threshold - 1);
      return detail::thresholdperiod(x * y, _threshold, _period);
    }

    //! \complexity
    //! Constant.
    Scalar plus_no_checks(Scalar x, Scalar y) const noexcept {
      LIBSEMIGROUPS_ASSERT(x >= 0 && x <= _period + _threshold - 1);
      LIBSEMIGROUPS_ASSERT(y >= 0 && y <= _period + _threshold - 1);
      return detail::thresholdperiod(x + y, _threshold, _period);
    }

    //! \complexity
    //! Constant.
    Scalar threshold() const noexcept {
      return _threshold;
    }

    //! \complexity
    //! Constant.
    Scalar period() const noexcept {
      return _period;
    }
 
   private:
    Scalar _period;
    Scalar _threshold;
  };

  template <typename Scalar>
  using DynamicNTPMatWithSemiring = DynamicMatrix<NTPSemiring<Scalar>, Scalar>;

  template <size_t T, size_t P, typename Scalar>
  using DynamicNTPMatWithoutSemiring = DynamicMatrix<NTPPlus<T, P, Scalar>,
                                                     NTPProd<T, P, Scalar>,
                                                     IntegerZero<Scalar>,
                                                     IntegerOne<Scalar>,
                                                     Scalar>;

  template <size_t T, size_t P, size_t R, size_t C, typename Scalar>
  using StaticNTPMat = StaticMatrix<NTPPlus<T, P, Scalar>,
                                    NTPProd<T, P, Scalar>,
                                    IntegerZero<Scalar>,
                                    IntegerOne<Scalar>,
                                    R,
                                    C,
                                    Scalar>;

  template <size_t T        = 0,
            size_t P        = 0,
            size_t R        = 0,
            size_t C        = R,
            typename Scalar = size_t>
  using NTPMat = std::conditional_t<
      R == 0 || C == 0,
      std::conditional_t<T == 0 && P == 0,
                         DynamicNTPMatWithSemiring<Scalar>,
                         DynamicNTPMatWithoutSemiring<T, P, Scalar>>,
      StaticNTPMat<T, P, R, C, Scalar>>;
 
  namespace detail {
    template <typename Scalar>
    struct IsNTPMatHelper<DynamicNTPMatWithSemiring<Scalar>> : std::true_type {
      static constexpr Scalar threshold = UNDEFINED;
      static constexpr Scalar period    = UNDEFINED;
    };
 
    template <size_t T, size_t P, typename Scalar>
    struct IsNTPMatHelper<DynamicNTPMatWithoutSemiring<T, P, Scalar>>
        : std::true_type {
      static constexpr Scalar threshold = T;
      static constexpr Scalar period    = P;
    };
 
    template <size_t T, size_t P, size_t R, size_t C, typename Scalar>
    struct IsNTPMatHelper<StaticNTPMat<T, P, R, C, Scalar>> : std::true_type {
      static constexpr Scalar threshold = T;
      static constexpr Scalar period    = P;
    };
 
    template <typename T>
    struct IsTruncMatHelper<T, std::enable_if_t<IsNTPMat<T>>> : std::true_type {
      static constexpr typename T::scalar_type threshold
          = IsNTPMatHelper<T>::threshold;
      static constexpr typename T::scalar_type period
          = IsNTPMatHelper<T>::period;
    };
  }  // namespace detail

  // Projective max-plus matrices
 
  namespace detail {
    template <typename T>
    class ProjMaxPlusMat : MatrixPolymorphicBase {
     public:
      using scalar_type            = typename T::scalar_type;
      using scalar_reference       = typename T::scalar_reference;
      using scalar_const_reference = typename T::scalar_const_reference;
      using semiring_type          = void;
 
      using container_type = typename T::container_type;
      using iterator       = typename T::iterator;
      using const_iterator = typename T::const_iterator;
 
      using underlying_matrix_type = T;
 
      using RowView = typename T::RowView;
 
      // Note that Rows are never normalised, and that's why we use the
      // underlying matrix Row type and not 1 x n ProjMaxPlusMat's instead
      // (since these will be normalised according to their entries, and
      // this might not correspond to the normalised entries of the matrix).
      using Row = typename T::Row;
 
      scalar_type scalar_one() const noexcept {
        return _underlying_mat.scalar_one();
      }
 
      scalar_type scalar_zero() const noexcept {
        return _underlying_mat.scalar_zero();
      }

      // ProjMaxPlusMat - Constructors + destructor - public
 
      ProjMaxPlusMat() : _is_normalized(false), _underlying_mat() {}
      ProjMaxPlusMat(ProjMaxPlusMat const&)            = default;
      ProjMaxPlusMat(ProjMaxPlusMat&&)                 = default;
      ProjMaxPlusMat& operator=(ProjMaxPlusMat const&) = default;
      ProjMaxPlusMat& operator=(ProjMaxPlusMat&&)      = default;
 
      ProjMaxPlusMat(size_t r, size_t c)
          : _is_normalized(false), _underlying_mat(r, c) {}
 
      // TODO(1) other missing constructors
      ProjMaxPlusMat(
          typename underlying_matrix_type::semiring_type const* semiring,
          size_t                                                r,
          size_t                                                c)
          : _is_normalized(false), _underlying_mat(semiring, r, c) {}
 
      explicit ProjMaxPlusMat(std::vector<std::vector<scalar_type>> const& m)
          : _is_normalized(false), _underlying_mat(m) {
        normalize();
      }
 
      ProjMaxPlusMat(
          std::initializer_list<std::initializer_list<scalar_type>> const& m)
          : ProjMaxPlusMat(
              std::vector<std::vector<scalar_type>>(m.begin(), m.end())) {}
 
      ~ProjMaxPlusMat() = default;
 
      ProjMaxPlusMat one() const {
        return ProjMaxPlusMat(_underlying_mat.one());
      }
 
      static ProjMaxPlusMat one(size_t n) {
        return ProjMaxPlusMat(T::one(n));
      }

      // Comparison operators
 
      bool operator==(ProjMaxPlusMat const& that) const {
        normalize();
        that.normalize();
        return _underlying_mat == that._underlying_mat;
      }
 
      bool operator!=(ProjMaxPlusMat const& that) const {
        return !(_underlying_mat == that._underlying_mat);
      }
 
      bool operator<(ProjMaxPlusMat const& that) const {
        normalize();
        that.normalize();
        return _underlying_mat < that._underlying_mat;
      }
 
      bool operator>(ProjMaxPlusMat const& that) const {
        return that < *this;
      }
 
      template <typename Thing>
      bool operator>=(Thing const& that) const {
        static_assert(IsMatrix<Thing> || std::is_same_v<Thing, RowView>);
        return that < *this || that == *this;
      }
 
      // not noexcept because operator< isn't
      template <typename Thing>
      bool operator<=(Thing const& that) const {
        static_assert(IsMatrix<Thing> || std::is_same_v<Thing, RowView>);
        return *this < that || that == *this;
      }

      // Attributes
 
      scalar_reference operator()(size_t r, size_t c) {
        // to ensure the returned value is normalised
        normalize();
        // to ensure that the matrix is renormalised if the returned scalar is
        // assigned.
        _is_normalized = false;
        return _underlying_mat(r, c);
      }
 
      scalar_reference at(size_t r, size_t c) {
        matrix::throw_if_bad_coords(*this, r, c);
        return this->operator()(r, c);
      }
 
      scalar_const_reference operator()(size_t r, size_t c) const {
        normalize();
        return _underlying_mat(r, c);
      }
 
      scalar_const_reference at(size_t r, size_t c) const {
        matrix::throw_if_bad_coords(*this, r, c);
        return this->operator()(r, c);
      }
 
      size_t number_of_rows() const noexcept {
        return _underlying_mat.number_of_rows();
      }
 
      size_t number_of_cols() const noexcept {
        return _underlying_mat.number_of_cols();
      }
 
      size_t hash_value() const {
        normalize();
        return Hash<T>()(_underlying_mat);
      }

      // Arithmetic operators - in-place
 
      void product_inplace_no_checks(ProjMaxPlusMat const& A,
                                     ProjMaxPlusMat const& B) {
        _underlying_mat.product_inplace_no_checks(A._underlying_mat,
                                                  B._underlying_mat);
        normalize(true);  // force normalize
      }
 
      void product_inplace(ProjMaxPlusMat const& A, ProjMaxPlusMat const& B) {
        _underlying_mat.product_inplace(A._underlying_mat, B._underlying_mat);
        normalize(true);  // force normalize
      }
 
      void plus_inplace_no_checks(ProjMaxPlusMat const& that) {
        _underlying_mat.plus_inplace_no_checks(that._underlying_mat);
        normalize(true);  // force normalize
      }
 
      void operator+=(ProjMaxPlusMat const& that) {
        _underlying_mat += that._underlying_mat;
        normalize(true);  // force normalize
      }
 
      void product_inplace_no_checks(ProjMaxPlusMat const& that) {
        _underlying_mat.product_inplace_no_checks(that._underlying_mat);
        normalize(true);  // force normalize
      }
 
      void operator*=(scalar_type a) {
        _underlying_mat *= a;
        normalize(true);  // force normalize
      }
 
      void operator+=(scalar_type a) {
        _underlying_mat += a;
        normalize(true);  // force normalize
      }
 
      void operator+=(RowView a) {
        _underlying_mat += a;
        normalize(true);  // force normalize
      }
 
      ProjMaxPlusMat operator*(scalar_type a) const {
        ProjMaxPlusMat result(*this);
        result *= a;
        return result;
      }
 
      ProjMaxPlusMat operator+(scalar_type a) const {
        ProjMaxPlusMat result(*this);
        result += a;
        return result;
      }

      // Arithmetic operators - not in-place
 
      ProjMaxPlusMat plus_no_checks(ProjMaxPlusMat const& that) const {
        return ProjMaxPlusMat(
            plus_no_checks(_underlying_mat, that._underlying_mat));
      }
 
      ProjMaxPlusMat operator+(ProjMaxPlusMat const& that) const {
        return ProjMaxPlusMat(_underlying_mat + that._underlying_mat);
      }
 
      ProjMaxPlusMat product_no_checks(ProjMaxPlusMat const& that) const {
        return ProjMaxPlusMat(
            product_no_checks(_underlying_mat, that._underlying_mat));
      }
 
      ProjMaxPlusMat operator*(ProjMaxPlusMat const& that) const {
        return ProjMaxPlusMat(_underlying_mat * that._underlying_mat);
      }

      // Iterators
 
      // The following should probably be commented out because I can't
      // currently think how to ensure that the matrix is normalised if it's
      // changed this way.
 
      iterator begin() noexcept {
        // to ensure the returned value is normalised
        normalize();
        // to ensure that the matrix is renormalised if the returned scalar is
        // assigned.
        _is_normalized = false;
        return _underlying_mat.begin();
      }
 
      iterator end() noexcept {
        // to ensure the returned value is normalised
        normalize();
        // to ensure that the matrix is renormalised if the returned scalar is
        // assigned.
        _is_normalized = false;
        return _underlying_mat.end();
      }
 
      const_iterator begin() const noexcept {
        normalize();
        return _underlying_mat.begin();
      }
 
      const_iterator end() const noexcept {
        normalize();
        return _underlying_mat.end();
      }
 
      const_iterator cbegin() const noexcept {
        normalize();
        return _underlying_mat.cbegin();
      }
 
      const_iterator cend() const noexcept {
        normalize();
        return _underlying_mat.cend();
      }

      // Modifiers
 
      void swap(ProjMaxPlusMat& that) noexcept {
        std::swap(_underlying_mat, that._underlying_mat);
      }
 
      void transpose() noexcept {
        _underlying_mat.transpose();
      }
 
      void transpose_no_checks() noexcept {
        _underlying_mat.transpose_no_checks();
      }

      // Rows
 
      RowView row(size_t i) const {
        normalize();
        return _underlying_mat.row(i);
      }
 
      template <typename C>
      void rows(C& x) const {
        normalize();
        return _underlying_mat.rows(x);
      }

      // Friend functions
 
      friend std::ostream& operator<<(std::ostream&         os,
                                      ProjMaxPlusMat const& x) {
        x.normalize();
        os << detail::to_string(x._underlying_mat);
        return os;
      }
 
      T const& underlying_matrix() const noexcept {
        normalize();
        return _underlying_mat;
      }
 
     private:
      explicit ProjMaxPlusMat(T&& mat)
          : _is_normalized(false), _underlying_mat(std::move(mat)) {
        normalize();
      }
 
      void normalize(bool force = false) const;
 
      mutable bool _is_normalized;
      mutable T    _underlying_mat;
    };
  }  // namespace detail


  template <size_t R, size_t C, typename Scalar>
  using StaticProjMaxPlusMat
      = detail::ProjMaxPlusMat<StaticMaxPlusMat<R, C, Scalar>>;

  template <typename Scalar>
  using DynamicProjMaxPlusMat
      = detail::ProjMaxPlusMat<DynamicMaxPlusMat<Scalar>>;

  template <size_t R = 0, size_t C = R, typename Scalar = int>
  using ProjMaxPlusMat = std::conditional_t<R == 0 || C == 0,
                                            DynamicProjMaxPlusMat<Scalar>,
                                            StaticProjMaxPlusMat<R, C, Scalar>>;
 
  namespace detail {
    template <size_t R, size_t C, typename Scalar>
    struct IsProjMaxPlusMatHelper<StaticProjMaxPlusMat<R, C, Scalar>>
        : std::true_type {};
 
    template <typename Scalar>
    struct IsProjMaxPlusMatHelper<DynamicProjMaxPlusMat<Scalar>>
        : std::true_type {};
  }  // namespace detail

  //! \no_libsemigroups_except
  //!
  //! \warning This function does not detect overflows of `Mat::scalar_type`.
  template <typename Mat>
  auto operator+(typename Mat::scalar_type a, Mat const& x)
      -> std::enable_if_t<IsMatrix<Mat>, Mat> {
    return x + a;
  }

  //! \no_libsemigroups_except
  //!
  //! \warning This function does not detect overflows of `Mat::scalar_type`.
  template <typename Mat>
  auto operator*(typename Mat::scalar_type a, Mat const& x)
      -> std::enable_if_t<IsMatrix<Mat>, Mat> {
    return x * a;
  }
}  // namespace libsemigroups
 
namespace std {
  template <size_t N,
            typename Mat,
            std::enable_if_t<libsemigroups::IsMatrix<Mat>>>
  inline void swap(Mat& x, Mat& y) noexcept {
    x.swap(y);
  }
}  // namespace std
 
#include "matrix-class.tpp"
#endif  // LIBSEMIGROUPS_MATRIX_CLASS_HPP_