// Diggins Public Domain Post #5
// Contributed to Public Domain by Christopher Diggins
// This is an efficient matrix class implementation for when rows and
columns are known at compile-time
#include <algorithm>
#include <iterator>
#include <cassert>
//#include "static_assert.hpp"
template<class Value_T, int Size_N, int Step_N>
class kslice_iter
{
public:
// public typedefs
typedef Value_T value_type;
typedef Value_T& reference;
typedef const Value_T& const_reference;
typedef std::size_t size_type;
typedef std::ptrdiff_t difference_type;
typedef kslice_iter self;
// constructors
explicit kslice_iter(Value_T* x) : m(x) {
}
// operators
self& operator++() {
m += step();
return *this;
}
self operator++(int) {
self tmp = *this;
m += step();
return tmp;
}
Value_T& operator[](size_type n) {
assert(n < size());
return m[n * Step_N];
}
const Value_T& operator[](size_type n) const {
assert(n < N);
return m[n * Step_N];
}
static int size() {
return Size_N;
}
static int step() {
return Step_N;
}
Value_T& operator*() {
return *m;
}
friend difference_type operator-(const self& x, const self& y) {
return (y.m - x.m) / step();
}
friend bool operator==(const self& x, const self& y) {
return x.m == y.m;
}
friend bool operator!=(const self& x, const self& y) {
return x.m != y.m;
}
friend bool operator<(const self& x, const self& y) {
return x.m < y.m;
}
private:
// default constructor hidden from view
kslice_iter() { };
Value_T* m;
};
template<int N, class T, class RowIter_T, class ColIter_T>
T dot_product(RowIter_T x, ColIter_T y)
{
T ret(0);
for (int i=N; i > 0; --i) {
ret += *x++ * *y++;
}
return ret;
}
template<class T, class M1, class M2, class M3>
kmatrix_multiply(const M1& m1, const M2& m2, M3& result)
{
//STATIC_ASSERT(M1::cols == M2::rows);
//STATIC_ASSERT(M3::rows == M1::rows);
//STATIC_ASSERT(M3::cols == M2::cols);
for (int i=0; i < M1::rows; ++i) {
for (int j=0; j < M2::cols; ++j) {
result[i][j] = dot_product<M1::cols, T>(m1.row(i), m2.column(j));
}
}
}
template<class Value_T, int Rows_N, int Cols_N>
class kmatrix
{
public:
// public typedefs
typedef Value_T value_type;
typedef kmatrix self;
typedef Value_T* iterator;
typedef const Value_T* const_iterator;
typedef kslice_iter<Value_T, Cols_N, 1> row_type;
typedef kslice_iter<Value_T, Rows_N, Cols_N> col_type;
// public constants
static const int rows = Rows_N;
static const int cols = Cols_N;
// constructors
kmatrix() {
std::fill(begin(), end(), Value_T(0));
}
kmatrix(Value_T& x) {
std::fill(begin(), end(), x);
}
kmatrix(const self& x) {
operator=(x);
}
// public functions
self& operator=(const self& x) {
std::copy(x.begin(), x.end(), begin());
return *this;
}
self& operator=(Value_T& x) {
std::fill(begin(), end(), x);
return *this;
}
row_type row(int n) const {
assert(n < rows);
return row_type(m + (n * Cols_N));
}
col_type column(int n) const {
assert(n < cols);
return col_type(m + n);
}
row_type operator[](int n) {
return row(n);
}
const row_type operator[](int n) const {
return row(n);
}
iterator begin() {
return m;
}
iterator end() {
return m + size();
}
const_iterator begin() const {
return m;
}
const_iterator end() const {
return m + size();
}
static int size() {
return Rows_N * Cols_N;
}
// operators
self& operator+=(const self& x) {
for (int i=0; i < size(); ++i) {
m[i] += x.m[i];
}
return *this;
}
self& operator-=(const self& x) {
for (int i=0; i < size(); ++i) {
m[i] -= x.m[i];
}
return *this;
}
self& operator+=(value_type x) {
for (int i=0; i < size(); ++i) {
m[i] += x;
}
return *this;
}
self& operator-=(value_type x) {
for (int i=0; i < size(); ++i) {
m[i] -= x;
}
return *this;
}
self& operator*=(value_type x) {
for (int i=0; i < size(); ++i) {
m[i] *= x;
}
return *this;
}
self& operator/=(value_type x) {
for (int i=0; i < size(); ++i) {
m[i] /= x;
}
return *this;
}
self& operator=(value_type x) {
std::fill(begin(), end(), x);
return *this;
}
self operator-() {
self x;
for (int i=0; i < size(); ++i) {
x.m[i] = -m[i];
}
return x;
}
// friends
friend self operator+(const self& x, const self& y) {
return self(x) += y;
}
friend self operator-(const self& x, const self& y) {
return self(x) -= y;
}
friend self operator+(const self& x, value_type y) {
return self(x) += y;
}
friend self operator-(const self& x, value_type y) {
return self(x) -= y;
}
friend self operator*(const self& x, value_type y) {
return self(x) *= y;
}
friend self operator/(const self& x, value_type y) {
return self(x) /= y;
}
template<int Cols2_N>
friend kmatrix<Value_T, Rows_N, Cols2_N>
operator*(const self& x, const kmatrix<Value_T, Cols_N, Cols2_N>& y) {
kmatrix<Value_T, Rows_N, Cols2_N> ret;
kmatrix_multiply(x, y, ret);
return ret;
}
private:
mutable Value_T m[Rows_N * Cols_N];
};
