#include <iostream>
#include <sstream>
#include <vector>

template <typename T> struct dense_traits;
template <typename T> struct dense_traits<const T> : dense_traits<T> {};
template <typename Parent, std::size_t M, std::size_t N> class Subview;
template <typename Derived> class Densebase;

template <typename Derived>
class Mathbase
{
public:
  using T = typename dense_traits<Derived>::T;

  Derived& derived() { return static_cast<Derived&>(*this); }
  const Derived& derived() const { return static_cast<const Derived&>(*this); }

  T& coeff(std::size_t row, std::size_t col) { return derived().coeff(row, col); }
  const T& coeff(std::size_t row, std::size_t col) const { return derived().coeff(row, col); }
  T& coeff(std::size_t index) { return derived().coeff(index); }
  const T& coeff(std::size_t index) const { return derived().coeff(index); }
};

template <typename Derived, std::size_t N>
struct det_helper
{
  static inline typename dense_traits<Derived>::T run(const Densebase<Derived>& dense)
  {
    typename dense_traits<Derived>::T det = 0;
    // Laplace expansion
    for (std::size_t i = 0; i < N; ++i)
      det += ((i & 1) ? -1 : 1)*dense.coeff(0,i)*det_helper::run(dense.sub(0,i));

    return det;
  }
};

template <typename Derived>
struct det_helper<Derived, 2>
{
  static inline typename dense_traits<Derived>::T run(const Densebase<Derived>& dense)
  {
    return dense.coeff(0,0)*dense.coeff(1,1) - dense.coeff(0,1)*dense.coeff(1,0);
  }
};

template <typename Derived>
class Densebase : public Mathbase<Derived>
{
public:
  using base = Mathbase<Derived>;
  using base::base;
  using base::derived;
  using base::coeff;
  using typename base::T;

  using Submat = Subview<Derived, dense_traits<Derived>::M-1, dense_traits<Derived>::N-1>;
  using ConstSubmat = const Subview<const Derived, dense_traits<Derived>::M-1, dense_traits<Derived>::N-1>;

  constexpr std::size_t rows() const { return derived().rows(); }
  constexpr std::size_t cols() const { return derived().cols(); }

  constexpr std::size_t diagonalsize() const { return std::min(rows(), cols()); }

  T determinant() const
  {
    return det_helper<Derived, dense_traits<Derived>::M>::run(derived());
  }

  Submat sub(std::size_t row, std::size_t col) { return Submat(derived(), row, col); }
  ConstSubmat sub(std::size_t row, std::size_t col) const { return ConstSubmat(derived(), row, col); }
};

template <typename T, std::size_t M, std::size_t N>
class Matrix : public Densebase<Matrix<T,M,N>>
{
public:
  Matrix() : data_(M*N) {}
  Matrix(std::initializer_list<T> vals) : data_{vals} {}
  virtual ~Matrix() {}

  T& coeff(std::size_t row, std::size_t col) { return data_[row*N + col]; }
  T& coeff(std::size_t index) { return data_[index]; }

  const T& coeff(std::size_t row, std::size_t col) const { return data_[row*N + col]; }
  const T& coeff(std::size_t index) const { return data_[index]; }

  constexpr std::size_t rows() const { return M; }
  constexpr std::size_t cols() const { return N; }

private:
  std::vector<T> data_;
};

template <typename Parent, std::size_t M, std::size_t N>
class Subview : public Densebase<Parent>
{
public:
  using T = typename dense_traits<Parent>::T;

  Subview(Parent& parent, std::size_t iM, std::size_t iN) : parent_{parent}, iM_{iM}, iN_{iN} {}
  T& coeff(std::size_t row, std::size_t col) { return parent_.coeff(row + (row >= iM_), col + (col >= iN_)); }
  T& coeff(std::size_t index) { return parent_.coeff(index + (index >= iM_), index + (index >= iN_)); }

  const T& coeff(std::size_t row, std::size_t col) const { return parent_.coeff(row + (row >= iM_), col + (col >= iN_)); }
  const T& coeff(std::size_t index) const { return parent_.coeff(index + (index >= iM_), index + (index >= iN_)); }

  constexpr std::size_t rows() const { return M; }
  constexpr std::size_t cols() const { return N; }

private:
  Parent& parent_;
  std::size_t iM_, iN_;
};

template <typename T_, std::size_t M_, std::size_t N_>
struct dense_traits< Matrix<T_, M_, N_> >
{
  using T = T_;
  enum { M = M_, N = N_ };
};

template <typename Derived>
inline std::ostream& operator<<(std::ostream& stream, const Densebase<Derived>& dense)
{
  unsigned int width = 0;
  for (std::size_t i = 0; i < dense.cols(); ++i)
    for (std::size_t j = 0; j < dense.rows(); ++j)
    {
      std::stringstream w; w.precision(3); w << std::fixed << dense.coeff(j,i);
      width = std::max<unsigned int>(width, w.tellp());
    }

  stream.precision(3);
  stream.setf(std::ios::fixed);
  for (std::size_t i = 0; i < dense.rows(); ++i)
  {
    stream << std::endl;
    stream.width(width); stream << dense.coeff(i,0);
    for (std::size_t j = 1; j < dense.cols(); ++j)
    {
      stream << ' ';
      stream.width(width); stream << dense.coeff(i,j);
    }
  }

  return stream;
}

int main(int argc, char* argv[])
{

  Matrix<float, 3, 3> m { 2, 1, 2,
                          2, 5, 4,
                          5, 6, 1};

  auto t = m.sub(0,0);
  std::cout << m.determinant();
  t.coeff(0,0) = 10;
  std::cout << m;
  std::cout << t.coeff(1,0);
  std::cout << m.sub(0,0);

  return EXIT_SUCCESS;
}