n-dimensional, generic, compile-time matrix and vector, reusing the same code

main.cpp

#include <iostream>
#include <cassert>
#include <cstring>
#include "nmath.h"

using namespace std;

typedef matrix_nxnt<2, float> matrix2x2;
typedef vector_nt<3, float> vector3;
typedef vector_nt<4, float> vector4;
typedef matrix_nxnt<3, float> matrix3x3;

int main()
{
    vector4 a = vector4::coord(1, 2, 3, 4);
    vector4 b = vector4::coord(2, 4, 6, 8);

    vector4 c = a + b;
    a += b;

    assert(a == c);
    assert(a == vector4::coord(3, 6, 9, 12));

    matrix2x2 ma = matrix2x2::rows(3, 5, 7, 9);
    matrix2x2 mb = matrix2x2::rows(1, 1, 1, 1);

    matrix2x2 mc = ma + mb;
    ma += mb;

    assert(ma == mc);
    assert(ma == matrix2x2::rows(4, 6, 8, 10));
    assert(ma.determinant() == -8);

    matrix3x3 m3 = matrix3x3::rows(
        1, 2, 3,
        5, 1, 1,
        5, 6, 7
    );

    assert(m3.determinant() == 16);

    assert(m3.row_vector(2) == vector3::coord(5, 6, 7));

    m3.row_vector(0) = vector3::coord(0, 0, 0);
    assert(m3.row_vector(0) == vector3::zero());

    vector3 v3ar[] = {
        vector3::coord(0, 0, 0),
        vector3::coord(5, 1, 1),
        vector3::coord(5, 6, 7),
    };

    float far[] = {
        0, 0, 0,
        5, 1, 1,
        5, 6, 7
    };

    assert(sizeof(m3) == sizeof(far));
    assert(sizeof(v3ar) == sizeof(far));
    assert(memcmp(&m3, far, sizeof(far)) == 0);
    assert(memcmp(v3ar, far, sizeof(far)) == 0);

    return 0;
}

nmath.h

#pragma once

template <typename _type, typename _this_type>
struct vector2_base
{
    typedef _type value_type;

    static _this_type coord(value_type x, value_type y)
    {
        _this_type ret;

        ret.x = x;
        ret.y = y;

        return ret;
    }

    _type x, y;
};

template <typename _type, typename _this_type>
struct vector3_base
{
    typedef _type value_type;

    static _this_type coord(value_type x, value_type y, value_type z)
    {
        _this_type ret;

        ret.x = x;
        ret.y = y;
        ret.z = z;

        return ret;
    }

    _type x, y, z;
};

template <typename _type, typename _this_type>
struct vector4_base
{
    typedef _type value_type;

    static _this_type coord(value_type x, value_type y, value_type z, value_type w)
    {
        _this_type ret;

        ret.x = x;
        ret.y = y;
        ret.z = z;
        ret.w = w;

        return ret;
    }

    _type x, y, z, w;
};

template <typename _type, typename _this_type>
struct matrix2x2_base
{
    typedef _type value_type;

    static _this_type rows(value_type rc_00, value_type rc_01, value_type rc_10, value_type rc_11)
    {
        _this_type ret;

        ret._00 = rc_00; ret._01 = rc_01;
        ret._10 = rc_10; ret._11 = rc_11;

        return ret;
    }

    _type _00, _01;
    _type _10, _11;
};

template <typename _type, typename _this_type>
struct matrix3x3_base
{
    typedef _type value_type;

    static _this_type rows(value_type rc_00, value_type rc_01, value_type rc_02,
                           value_type rc_10, value_type rc_11, value_type rc_12,
                           value_type rc_20, value_type rc_21, value_type rc_22)
    {
        _this_type ret;

        ret._00 = rc_00; ret._01 = rc_01; ret._02 = rc_02;
        ret._10 = rc_10; ret._11 = rc_11; ret._12 = rc_12;
        ret._20 = rc_20; ret._21 = rc_21; ret._22 = rc_22;

        return ret;
    }

    _type _00, _01, _02;
    _type _10, _11, _12;
    _type _20, _21, _22;
};



template <size_t n, typename _type, typename _this_type>
struct order_traits
{
};

template <typename _type, typename _this_type>
struct order_traits<2, _type, _this_type>
{
    typedef vector2_base<_type, _this_type> vector;
    typedef matrix2x2_base<_type, _this_type> matrix;
};

template <typename _type, typename _this_type>
struct order_traits<3, _type, _this_type>
{
    typedef vector3_base<_type, _this_type> vector;
    typedef matrix3x3_base<_type, _this_type> matrix;
};


template <typename _type, typename _this_type>
struct order_traits<4, _type, _this_type>
{
    typedef vector4_base<_type, _this_type> vector;
};


template <size_t _count, typename _this_type, typename _parent_type>
class ntuple : public _parent_type
{
public:
    typedef typename _parent_type::value_type value_type;
    typedef size_t size_type;
    static const size_type value_count = _count;

    static _this_type zero()
    {
        _this_type ret;

        //could use memset, but it's probably slower...
        //for should be optimized at compile time.
        //too lazy to test
        for(size_type i=0; i<value_count; ++i)
        {
            ret.at(i) = value_type(0);
        }

        return ret;
    }

    static _this_type& attach_to_ptr(value_type* ptr)
    {
        return *reinterpret_cast<_this_type*>(ptr);
    }

    static const _this_type& attach_to_ptr(const value_type* ptr)
    {
        return *reinterpret_cast<const _this_type*>(ptr);
    }

    value_type* as_array()
    {
        return reinterpret_cast<value_type*>(this);
    }

    const value_type* as_array() const
    {
        return reinterpret_cast<const value_type*>(this);
    }

    value_type& at(size_type i)
    {
        return as_array()[i];
    }

    const value_type& at(size_type i) const
    {
        return as_array()[i];
    }

    value_type& operator[](size_type i)
    {
        return at(i);
    }

    const value_type& operator[](size_type i) const
    {
        return at(i);
    }

    _this_type& operator+=(const _this_type& b)
    {
        for(size_type i=0; i<value_count; ++i)
        {
            at(i) += b.at(i);
        }

        return as_this_type();
    }

protected:
    _this_type& as_this_type()
    {
        return *reinterpret_cast<_this_type*>(this);
    }

    const _this_type& as_this_type() const
    {
        return *reinterpret_cast<const _this_type*>(this);
    }
};

template <size_t _count, typename _this_type, typename _parent_type>
bool operator==(const ntuple<_count, _this_type, _parent_type>& a,
                const ntuple<_count, _this_type, _parent_type>& b)
{
    for(size_t i=0; i<_this_type::value_count; ++i)
    {
        if(a.at(i) != b.at(i))
            return false;
    }

    return true;
}

template <size_t _count, typename _this_type, typename _parent_type>
bool operator!=(const ntuple<_count, _this_type, _parent_type>& a,
                const ntuple<_count, _this_type, _parent_type>& b)
{
    return !operator==(a, b);
}


template <size_t _count, typename _this_type, typename _parent_type>
_this_type operator+(const ntuple<_count, _this_type, _parent_type>& a,
                     const ntuple<_count, _this_type, _parent_type>& b)
{
    _this_type ret;

    for(size_t i=0; i<_this_type::value_count; ++i)
        ret.at(i) = a.at(i) + b.at(i);

    return ret;
}

template <size_t _dim, typename _type>
class vector_nt : public ntuple<_dim, vector_nt<_dim, _type>, typename order_traits<_dim, _type, vector_nt<_dim, _type> >::vector>
{
    typedef ntuple<_dim, vector_nt<_dim, _type>, typename order_traits<_dim, _type, vector_nt<_dim, _type> >::vector> super;
public:
    typedef typename super::size_type size_type;
    using super::value_count;


    static const size_type dimension = value_count;
};

template <size_t _dim, typename _type>
std::ostream& operator<<(std::ostream& o, const vector_nt<_dim, _type>& vec)
{
    typedef vector_nt<_dim, _type> vector;

    o << '(';
    for(typename vector::size_type i=0; i<vector::dimension; ++i)
    {
        if(i)
        {
            o << ", ";
        }

        o << vec.at(i);
    }
    o << ')';

    return o;
}

template <size_t _n, typename _type>
class matrix_nxnt : public ntuple<_n*_n, matrix_nxnt<_n, _type>, typename order_traits<_n, _type, matrix_nxnt<_n, _type> >::matrix>
{
    typedef ntuple<_n*_n, matrix_nxnt<_n, _type>, typename order_traits<_n, _type, matrix_nxnt<_n, _type> >::matrix> super;
public:
    typedef typename super::size_type size_type;
    typedef typename super::value_type value_type;

    static const size_type order = _n;

    value_type* row(size_type i)
    {
        return this->as_array() + i*order;
    }

    const value_type* row(size_type i) const
    {
        return this->as_array() + i*order;
    }

    value_type& operator()(size_type row, size_type column)
    {
        return this->row(row)[column];
    }

    const value_type& operator()(size_type row, size_type column) const
    {
        return this->row(row)[column];
    }

    value_type adjunct(size_type row, size_type col) const
    {
         return first_minor(row, col) * (((row+col)&1)?value_type(-1):value_type(1));
    }

    value_type first_minor(size_type row, size_type col) const
    {
        matrix_nxnt<order-1, value_type> mnr;
        for(size_type r=0; r<row; ++r)
        {
            for(size_type c=0; c<col; ++c)
                mnr(r, c) = (*this)(r, c);

            for(size_type c=col+1; c<order; ++c)
                mnr(r, c-1) = (*this)(r, c);
        }

        for(size_type r=row+1; r<order; ++r)
        {
            for(size_type c=0; c<col; ++c)
                mnr(r-1, c) = (*this)(r, c);

            for(size_type c=col+1; c<order; ++c)
                mnr(r-1, c-1) = (*this)(r, c);
        }

        return mnr.determinant();
    }

    vector_nt<order, value_type>& row_vector(size_type i)
    {
        return vector_nt<order, value_type>::attach_to_ptr(row(i));
    }

    const vector_nt<order, value_type>& row_vector(size_type i) const
    {
        return vector_nt<order, value_type>::attach_to_ptr(row(i));
    }


    value_type determinant() const
    {
        return det_impl<order>::eval(*this);
    }

private:

    //_unused is a stupid gcc workaround
    //it doesn't allow fully specialized member types, or any specialization of methods
    //
    //update:
    //As http://tinyurl.com/bmaa9yh suggests, this is (unexplicably) fobidden by the standard.
    //One would wonder why this is so, considering that it works fine with msvc without the dummy type
    //and even as a specialized template method.
    template <size_type _d, typename _unused = int>
    struct det_impl
    {
        static value_type eval(const matrix_nxnt<_d, value_type>& m)
        {
            //laplace
            value_type det = 0;
            for(size_type i=0; i<_d; ++i)
            {
                det += m(0, i) * m.adjunct(0, i);
            }


            return det;
        }
    };

    template <typename _unused>
    struct det_impl<2, _unused>
    {
        static value_type eval(const matrix_nxnt<2, value_type>& m)
        {
            return m._00*m._11 - m._01*m._10;
        }
    };
};

template <size_t _n, typename _type>
std::ostream& operator<<(std::ostream& o, const matrix_nxnt<_n, _type>& m)
{
    typedef matrix_nxnt<_n, _type> matrix;

    for(typename matrix::size_type r=0; r<matrix::order; ++r)
    {
        if(r)
        {
            o << std::endl;
        }

        o << '(';
        for(typename matrix::size_type c=0; c<matrix::order; ++c)
        {
            if(c)
            {
                o << ", ";
            }

            o << m(r, c);
        }
        o << ')';
    }

    return o;
}