I am implementing a litte generic math library. What I have done is to write my generic matrix and vector class. I'm curious if I have it done right so far (implementation wise not totally mathematical correctness wise). Reason: I'm relatively new to template programming and I am not sure what is right, what could be done better and what should I try to avoid at all costs if I use templates in C++.
So far I managed to multiply different matrices of (obviously) different sizes while maintaining static allocation by using only std::array
, so no dynamic memory allocation.
One thing is needed to be mentioned: my Vectors
and Matrices
start count their elements with #1 not #0 like it is usually done in C++ with arrays. The reason for this is to be near as possible to 'proper' math, but I am open to get convinced otherwise.
First of all, my generic Vector
class is quite boring but needed for my matrix class, so I put it inside this post.
#include <array>
#include <cmath>
namespace jslmath
{
template<size_t Dimension,typename NumberType = double>
class Vector
{
public:
using Value = NumberType;
using Storage = std::array<NumberType, Dimension>;
private:
Storage mField;
public:
Vector(const Vector&) = default;
Vector(Vector&&) = default;
virtual ~Vector() = default;
template<typename ...Targs>
Vector(Targs... args): mField({args...}){}
Vector(const Storage& args) : mField(args){}
Vector& operator=(const Vector&) = default;
Vector& operator=(Vector&&) = default;
Value operator[](size_t index)const
{
return at(index);
}
Value at(size_t index)const
{
if (index <= Dimension && index != 0)
return mField[index - 1];
throw;//need improvment
}
Vector operator+(Vector vec) const{ return add(vec); }
Vector operator-(Vector vec) const{ return sub(vec); };
Vector operator/(Vector vec) const{ return div(vec); };
Vector operator*(Vector vec) const{ return mul(vec); };
Vector operator+(Value scalar) const{ return add(scalar); };
Vector operator-(Value scalar) const{ return sub(scalar); };
Vector operator/(Value scalar) const{ return div(scalar); };
Vector operator*(Value scalar) const{ return mul(scalar); };
Vector add(Vector vec)const
{
Storage temp;
for (auto i = 0; i < Dimension; i++)
temp[i] = mField[i] + vec[i + 1];
return {temp};
}
Vector sub(Vector vec)const
{
Storage temp;
for (auto i = 0; i < Dimension; i++)
temp[i] = mField[i] - vec[i + 1];
return { temp };
}
Vector div(Vector vec)const
{
Storage temp;
for (auto i = 0; i < Dimension; i++)
temp[i] = mField[i] / vec[i + 1];
return { temp };
}
Vector mul(Vector vec)const
{
Storage temp;
for (auto i = 0; i < Dimension; i++)
temp[i] = mField[i] * vec[i + 1];
return { temp };
}
Value dot(Vector vec) const
{
Value tmp = 0;
for (auto i = 1; i <= Dimension; i++)
tmp = vec[i] + at(i);
return tmp;
}
Value magnitude() const
{
return std::sqrt(dot(*this));
}
Value magnitudeSq() const
{
return (dot(*this));
}
Value distance(Vector vec) const
{
auto tmp = vec - *this;
return tmp.magnitude();
}
void normalize()
{
*this = *this * (1.0 / magnitude());
}
double Angle(Vector vec)
{
return std::acos(dot(vec)/ std::sqrt(magnitudeSq()*vec.magnitudeSq()));
}
Vector add(Value scalar)const
{
Storage temp;
for (auto i = 0; i < Dimension; i++)
temp[i] = mField[i] + scalar;
return { temp };
}
Vector sub(Value scalar)const
{
Storage temp;
for (auto i = 0; i < Dimension; i++)
temp[i] = mField[i] - scalar;
return { temp };
}
Vector div(Value scalar)const
{
Storage temp;
for (auto i = 0; i < Dimension; i++)
temp[i] = mField[i] / scalar;
return { temp };
}
Vector mul(Value scalar)const
{
Storage temp;
for (auto i = 0; i < Dimension; i++)
temp[i] = mField[i] * scalar;
return { temp };
}
};
So from here on now the interesting part starts:
template<size_t N,size_t M,typename NumberType = double>
class Matrix
{
public:
using Value = NumberType;
using Storage = std::array<NumberType, N*M>;
using RowVec = Vector<M, NumberType>;
using ColVec = RowVec;
The first part of my class is just some lifetime saver and also makes the code more readable:
private:
Storage mGrid;
public:
Matrix(const Matrix&) = default;
Matrix(Matrix&&) = default;
template<typename ...Targs>
Matrix(Targs... args) : mGrid({ args... }) { }
Matrix(const Storage& args) : mGrid(args) { }
RowVec operator[](size_t index) const
{
return Row(index);
}
Value operator()(size_t Row, size_t Col) const
{
if(Row != 0 && Col != 0)
return mGrid[(Col-1) + M*(Row-1)];
throw;
}
constexpr size_t Height() { return N; };
constexpr size_t Width() { return M; }
Storage& data() { return mGrid; }
RowVec Row(size_t index)const
{
if (index <= N && index != 0)
{
typename RowVec::Storage temp;
for (auto i = 1; i <= M; i++)
temp[i - 1] = (*this)(index, i);
return RowVec(temp);
}
throw;
}
RowVec Col(size_t index)const
{
if (index <= N && index != 0)
{
typename RowVec::Storage temp;
for (auto i = 1; i <= N; i++)
temp[i - 1] = (*this)(i, index);
return RowVec(temp);
}
throw;
}
void transpose(){
Storage tmp;
for (auto i = 0; i < N * M; i++) {
int row = i / N;
int col = i % M;
tmp[i] = mGrid[M * col + row];
}
mGrid = tmp;
}
template< template<size_t, size_t, typename>class B, size_t I, size_t J, typename Type>
auto operator*(B<I, J, Type>& b)
{
return mul(b);
}
auto operator*(Value val)
{
return mul(val);
}
The most complicated part of my code so far. I have to calculate the size of my new matrix at compile time so that I am able to keep my goal to avoid dynamic allocation. The syntax of template templates is quite odd to me but it does work.
template< template<size_t,size_t,typename>class B, size_t I, size_t J, typename Type>
auto mul(B<I,J,Type>& b) -> decltype(Matrix<N, J, Value>{})
{
Matrix<N, J, Value> result;
for (auto i = 0; i < N; ++i)
for (auto j = 0; j < J; ++j)
{
for (int k = 0; k < I; ++k)
{
int _a = M * i + k;
int _b = J * k + j;
result.data()[J * i + j] += this->mGrid[_a] * b.data()[_b];
}
}
return result;
}
Matrix mul(Value b)
{
Storage result;
for (auto i = 0; i < N; ++i)
for (auto j = 0; j < M; ++j)
{
result[M*i + j] = mGrid[M*i + j] * b;
}
return { result };
}
};
}
And a small test application:
int main(int argc, char* argv[])
{
jslmath::Matrix<3, 3> Test3( 1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0 );
jslmath::Matrix<3, 1> Test4( 1.0, 2.0, 3.0 );
auto T = Test3 * Test4;
return T(1,1);
}
With the help of an online compiler I already found out that VS is not particularly great at compiling this code example compared to gcc and clang. What totally surprised me was that clang was so good at optimizing my code that in the and all that was left were a single return statement with value 14!
Clang did all the computation of multiplying two different sized matrices at compile time. I am flabbergasted about that. It is also the reason why I posted this code here. I want to know from you what I did well and what could be improved. I'm not really sure how I did that on the first try without even thinking of active optimization from my side.