mr-edd.co.uk :: horsing around with the C++ programming language

Filling matrices and vectors naturally

[4th November 2007]

The mathematics library I'm working on deals with vectors and matrices whose dimensions are given at compile time as template parameters. Ideally, for any given matrix or vector, we'd like a constructor that takes the corresponding number of elements in order to fill it directly. Unfortunately, we can't do this yet in C++. When initializer lists (Google Tech-talk video) are added to the next revision of the C++ standard, this problem will hopefully disappear, but I still need something in the mean time. This post is about my search for that something. I could have quite easily called this Why initializer lists are needed!

The problem

Here's a sketch of the problem and a handful of candidate solutions:

template<typename Element, unsigned N>
class vector
{
    public:
        // how to add a constructor that takes N elements?

        vector(...); // (1) varargs

        vector(const Element elements[N]); // (2) read from an array

        template<typename InputIter>
        vector(InputIter begin, InputIter end); // (3) read from an iterator range

        vector(const arg_container<Element> &args); // (4) use a custom type (we'll see how this works later)

    private:
        Element elements_[N];
};

I've labelled four candidate solutions. They each have their own problems, which we'll look at in turn.

varargs …

The first solution, (1), is to use varargs, which was originally designed for use in C in functions such as printf() and scanf().

While it can indeed provide a convenient way to specify an arbitrary number of arguments to a function or constructor, it simply isn't safe. Since no types are specified, no type checking can be done by the compiler. Furthermore, does the compiler call the copy-constructor of the objects passed as arguments? Or are they passed as references? I honestly don't know, and I don't care to know either. It's a dark corner of the language that was inherited from C for compatibility reasons and should be ignored in C++. It's hard to think of any good caveats[1]. Even if we use types without copy constructors, we still run in to problems:

typedef vector<float, 3> vec3f; // note the float!

vec3f u(12, -1.0, 3.14); // oops! doubles and an int

Here we passed in doubles and ints to the constructor, when we should have given floats. On platforms where sizeof(double) != sizeof(float) and/or sizeof(double) != sizeof(int), which are very common, this innocent little snippet leads to all kinds of disaster.

Furthermore, not only do we throw type checking out the window, but the compiler also has no idea how many arguments to expect. We know to pass N, but the compiler has no way of checking that N were actually given. The following code will compile just fine, but sends us to the dreaded plains of undefined behaviour as the constructor will access a third argument that wasn't really supplied.

typedef vector<double, 3> vec3;

vec3 u(-1.0, 3.14); // oops! 

So I abandoned the varargs idea pretty quickly[2].

Passing in an array

Solution number (2) is to have a constructor that takes an array of arguments. There are two problems with this. The first is that we might not always have an array of arguments available and we'd have to construct one just to pass to the constructor. The second problem, which is more serious, is that even though the signature of the constructor specifies an array of N elements, we can actually pass in any pointer-to-Element object. Again there's no way for the compiler to check that the array is of the correct size.

typedef vector<double, 3> vec3;

const double a[2] = {-1.0, 3.14};

vec3 u(a); // oops! 

The above will compile just fine but again, we have undefined behaviour as the constructor thinks we've passed in an array of 3 elements and so it will blindly access a[2]. So while we've eliminated the type-safety problem we saw in solution (1), we still aren't allowing the compiler to check the number of arguments given.

However, we can improve this situation a little by changing the constructor's signature to take a reference to an array:

template<typename Element, unsigned N>
class vector
{
    public:
        vector(const Element (&elements)[N]); // (2.1) read from a reference to an array
};

This now means that the compiler can enforce that an array with N elements is really given. Here's a quick example:

void give_me_an_array(const double (&a)[3])
{
}

int main()
{
    double a1[3] = { 1.0, 2.0, 3.0 };
    double a2[2] = { 1.0, 2.0 };

    give_me_an_array(a1); // fine
    give_me_an_array(a2); // compiler error

    return 0;
}

We still have to construct an array in order to call this constructor, but it is certainly much better than what we've seen so far.

Passing a pair of iterators

The third option is to pass in a pair of iterators. This means that we can use a pair of pointers (from an array), or elements of an STL container, or a number of other things. While this solution is certainly the most flexible seen so far we lose the compile-time bounds checking of the reference-to-an-array approach:

typedef vector<double, 2> vec2;

double a[3] = { 1.0, 2.0, 3.0 };
vec2 u(a, a + 3); // oops

However, we can claw something back with a run-time check inside the constructor. It's not ideal, but I would argue that not much is lost in this particular situation:

template<typename Element, unsigned N>
class vector
{
    public:
        template<typename InputIter>
        vector(InputIter begin, InputIter end)
        {
            assert(std::difference(begin, end) == N);
            std::copy(begin, end, elements_);
        }

    private:
        Element elements_[N];
};

But we still don't have the convenience of being able to pass the elements to the constructor directly e.g. vec3 u(1.0, 2.0, 3.0);

Creating a new object to allow a concise inline syntax

The 4th option is to invent a type that provides a convenient inline syntax that we can use to fill our vector:

typedef vector<double, 4> vec4;

// Create a 4D vector whose elements are 1, 2, 3 and 4
vec4 u( args<double>(1)(2)(3)(4) );

This is inspired by Andrei Alexandrescu's CUJ article, Inline containers for variable arguments. Here's a complete compilable implementation:

template<typename Element>
class arg_container
{
    public:
        typedef typename std::vector<Element>::size_type size_type;

        arg_container(Element e)
        {
            elements_.reserve(4); // vectors of size 2, 3 and 4 are most common for me
            elements_.push_back(e);
        }

        arg_container &operator() (Element e)
        {
            elements_.push_back(e);
            return *this;
        }

        size_type size() const
        {
            return elements_.size();
        }

        const Element *elements() const
        {
            return  elements_.empty() ?
                    static_cast<Element *>(0) :
                    &elements_.front();
        }

    private:
        std::vector<Element> elements_;
};

template<typename Element>
arg_container<Element> args(const Element &e)
{
    return arg_container<Element>(e);
}

template<typename Element, unsigned N>
class vector
{
    public:
        vector(const arg_container<Element> &ac)
        {
            typename arg_container<Element>::size_type sz = ac.size();
            assert(sz <= N);
            const Element *begin = ac.elements();
            std::copy(begin, begin + sz, elements_);
            std::fill(elements_ + sz, elements_ + N, Element(0));
        }

        Element operator[] (unsigned idx) const
        {
            assert(idx < N);
            return elements_[idx];
        }

    private:
        Element elements_[N];
};

int main()
{
    vector<double, 4> u( args<double>(1)(2)(3)(4) );

    for (unsigned i = 0; i != 4; ++i)
        std::cout << u[i] << '\\n';

    return 0;
}

While this effectively gives us a safe replacement for the varargs approach, it is perhaps a little ugly. Note that we have to explicitly specify the type double when we construct the vector, else the types for the arguments are deduced from the first element given, which would be int in my example code. We could arguably improve things a little by adapting the code to use an alternative syntax, by overloading the assignment and comma operators for arg_container:

typedef vector<double, 4> vec4;

// Create a 4D vector whose elements are 1, 2, 3 and 4
vec4 u( args<double>() = 1, 2, 3, 4 );

Even if this is better, we're doing some dynamic allocation here. vector construction is the kind of thing that's performed a lot inside tight inner loops, so ideally we'd like to avoid the overhead incurred by using std::vector<> underneath. The whole thing just feels a little bit heavy for the construction of a simple concrete vector class, anyway. Lastly, the bounds check is performed at run-time by an assert(), which is less desirable than having the check at compile-time.

It may be possible to use the stack instead of the heap and add an extra integral template argument to args_container<> in order to check that the number of elements given does not exceed a given value. However, I can't see a way of doing this right now.

Reaching a compromise

It's a shame that there is no perfect solution to this little problem[3]. So at some point a there has to be a compromise. Let's quickly summarise the properties of our ideal constructor(s):

  1. We should be able to construct a vector or matrix by passing the elements directly to the constructor
  2. Avoid dynamic allocation
  3. Concise and natural syntax
  4. Compile time bounds checking

None of the techniques we've seen so far has all of these features. However, I realised that in the vast majority of cases, I was wanting to create 2D, 3D or 4D vectors. This means that by providing constructors in my vector class that take 2, 3 and 4 elements, I could have my cake and eat it for small sized vectors:

template<typename Element, unsigned N>
class vector
{
    public:
        vector(Element x, Element y)
        {
            BOOST_STATIC_ASSERT(N == 2);
            elements_[0] = x; elements_[1] = y;
        } 

        vector(Element x, Element y, Element z)
        {
            BOOST_STATIC_ASSERT(N == 3);
            elements_[0] = x; elements_[1] = y; elements_[2] = z;
        }  

        vector(Element x, Element y, Element z, Element t)
        {
            BOOST_STATIC_ASSERT(N == 4);
            elements_[0] = x; elements_[1] = y; elements_[2] = z; elements_[3] = t;
        }  

    private:
        Element elements_[N];
};

Since my vector<> is a template class, it means that its member functions (including the constructors) are only instantiated if they're used. This means that the constructor taking two Elements should never get called for a 3D vector. But if someone does attempt to call it, it will fail to compile because of the static assertions.

typedef vector<double, 3> vec3;

vec3 u(1.0, 2.0, 3.0); // fine
vec3 u(1.0, 2.0); // compiler error

This still doesn't solve the problem when I want to construct vectors containing more than 4 elements of course, but it covers 90% of my use cases. For the remaining 10%, it is often the case that the elements are already sitting inside an array or container of some kind. We can use some of the techniques we're already seen to cover those cases.

But I was still left with the situation where I had to fill something like a 4×3 matrix a bunch of times in my test cases, for example. Without a nice way of doing this, I'm sure the tedium would have got to me! So what did I do? A bunch of existing libraries for linear algebra and mathematics (e.g. Blitz++) overload the assignment and comma operators so that you can fill vectors and matrices like this:

matrix3x3 m;
m = 1, 2, 3,
    4, 5, 6,
    7, 8, 9;

We can't pass elements to the constructor directly, but it does provide us with a convenient and natural way to fill the matrix. It also avoids dynamic allocation, which is a plus. But can we get compile-time bounds checking? As it turns out, the answer is yes, we can!

I felt a bit uncomfortable overloading the operator= to accept an element, so I went a slightly different route:

typedef matrix<double, 3, 3> mat3x3
mat3x3 m;
fill(m) = 1, 2, 3,
          4, 5, 6,
          7, 8, 9;

Here, the fill() function returns a type that has the appropriate operator= to begin assigning elements to the matrix. This keeps the matrix and vector interfaces free of any oddities. In fact, by using this approach I was able to make fill() work on native 1D and 2D arrays because 1D and 2D arrays fulfil the requirements to be models of the FixedVector and FixedMatrix concepts in my library:

double m[3][3];
fill(m) = 1, 2, 3,
          4, 5, 6,
          7, 8, 9;

Taking a look under the hood

So what's the code behind this and how is compile time bounds checking achieved?

The basic idea is that fill(m) returns an object that refers to the matrix m and is prepared to fill its first element when something is assigned to it. When this assignment happens, an object is returned whose comma operator is overloaded to fill the next element of m, and so on… By keeping track of the number of elements filled so far in a template parameter, we can enable compile time bounds checking. Here's the code straight out of my library:

namespace detail
{
    /*! This is used as a base class for comma_initialiser_hook. There are two specialisations for when Category
     *  is matrix_tag and vector_tag. In each specialisation, a static set() function is defined that sets the
     *  element at position Next to a particular value.
     */
    template<typename VecOrMat, typename Category, std::size_t Next>
    class element_setter;

    template<typename MutableFixedMatrix, std::size_t Next>
    class element_setter<MutableFixedMatrix, matrix_tag, Next>
    {
        private:
            static const std::size_t j = Next % matrix_width<MutableFixedMatrix>::value;
            static const std::size_t i = Next / matrix_width<MutableFixedMatrix>::value;

        protected:
            ~element_setter() { }

            static void set(MutableFixedMatrix &m, typename element_type<MutableFixedMatrix>::type e)
            {
                BOOST_STATIC_ASSERT(j < matrix_width<MutableFixedMatrix>::value);
                BOOST_STATIC_ASSERT(i < matrix_height<MutableFixedMatrix>::value);
                at(m, i, j) = e;
            }
    };

    template<typename MutableFixedVector, std::size_t Next>
    class element_setter<MutableFixedVector, vector_tag, Next>
    {
        protected:
            ~element_setter() { }

            static void set(MutableFixedVector &v, typename element_type<MutableFixedVector>::type e)
            {
                BOOST_STATIC_ASSERT(Next < vector_size<MutableFixedVector>::value);
                at(v, Next) = e;
            }
    };

    /*! This is the type returned from fill() in order to provide a means for filling matrices and vectors using a
     *  natural comma-based syntax.
     */
    template<typename MutableFixedThing, std::size_t Next>
    class comma_initializer_hook :
        private element_setter<MutableFixedThing, typename operand_category<MutableFixedThing>::type, Next>
    {
        private:
            typedef comma_initializer_hook<MutableFixedThing, Next + 1> next_hook;
            typedef typename element_type<MutableFixedThing>::type element_type;

        public:
            /*! Creates an object that is ready to fill the given vector or matrix */
            explicit comma_initializer_hook(MutableFixedThing *mft) : mft_(mft) { }

            /*! Assigns e to the first element of the thing that was passed in at construction and returns
             *  an object whose comma operator is overloaded to fill the element at position Next + 1.
             */
            next_hook operator= (element_type e)
            {
                BOOST_STATIC_ASSERT(Next == 0);
                set(*mft_, e);
                return next_hook(mft_);
            }

            /*! Assigns e to the Next'th element of the thing that was passed in at construction and returns
             *  an object whose comma operator is overloaded to fill the element at position Next + 1.
             */
            next_hook operator, (element_type e)
            {
                BOOST_STATIC_ASSERT(Next > 0);
                set(*mft_, e);
                return next_hook(mft_);
            }

        private:
            comma_initializer_hook &operator= (const comma_initializer_hook &); // remains undefined

            MutableFixedThing * const mft_;
    };

} // close namespace detail

template<typename MutableFixedThing>
detail::comma_initializer_hook<MutableFixedThing, 0>
fill(MutableFixedThing &mft)
{
    return detail::comma_initializer_hook<MutableFixedThing, 0>(&mft);
}

Here, I'm using a host of traits classes to get compile time information about matrix and vector types such as the dimensions (matrix_width<>, matrix_height<> and vector_size<>) and whether a given type is a matrix or a vector (operand_category<>). For more information on these traits, see this post.

But you can see that the overloaded assignment and comma operators return an object where the next element to be assigned is encoded as a template parameter, allowing the use of BOOST_STATIC_ASSERT() to do the bounds checking at compile time.

Conclusion

It's amazing how something that seems simple on the face of things actually turns out to be so fiddly. I'm hoping that with the introduction of initializer lists in the next C++ standard, we'll have a much cleaner way to initialise classes such as matrix and vector implementations.

I'd still like to be able to write something like:

typedef vector<double, 5> vec5;

vec5 u(args() = 1, 2, 3, 4, 5);

If anyone can see an efficient and type-safe way of doing so, please let me know! It might be possible to use Boost.Fusion to create a tuple, but I haven't got around to researching this just yet.

Footnotes
  1. The only valid use I've seen is in template meta-functions that check for convertibility between types []
  2. actually that's a lie! I didn't even consider it at the time because the associated problems are so well known. I only thought about it properly when writing this post! []
  3. at least none that I can see — if you know different please let me know! []

Comments

maxim1000

[03/06/2008 at 11:01:00]

Here is my try :)
The idea is that there is a chain of objects collecting elements and at the same time deducing some common type for them.
If I don’t mistake there was some implementation of extracting common type from two ones using ternary operation, but currently I don’t remember neither the way nor the place I saw this, so here dumb implementation is used.

template<typename type1,typename type2>
struct CommonType
{
};
template<typename type>
struct CommonType<type,type>
{
    typedef type Type;
};
template<>
struct CommonType<int,double>
{
    typedef double Type;
};
template<>
struct CommonType<double,int>
{
    typedef double Type;
};

template<typename commonType,unsigned currentDimension>
struct InlineContainer
{
    commonType Elements[currentDimension];

    template<typename newType>
    InlineContainer<typename CommonType<commonType,newType>::Type,currentDimension+1> operator<<(const newType &newElement)
    {
        InlineContainer<CommonType<commonType,newType>::Type,currentDimension+1> result;
        for(int c=0;c<currentDimension;++c)
            result.Elements[c]=Elements[c];
        result.Elements[currentDimension]=newElement;
        return result;
    }
};
struct ToStart
{
    template<typename newType>
    InlineContainer<newType,1> operator<<(const newType &newElement)
    {
        InlineContainer<newType,1> result;
        result.Elements[0]=newElement;
        return result;
    }
};

template<typename type,unsigned dimension>
class Vector
{
public:
    template<typename type2>
    Vector(const InlineContainer<type2,dimension> &container)
    {
        for(unsigned c=0;c<dimension;++c)
            Contents[c]=container.Elements[c];
    }
public:
    type Contents[dimension];
};

...
Vector<double,3> vector(ToStart()<<1<<2.5<<3);//fine
Vector<double,4> vector(ToStart()<<1<<2.5<<3);//error
Vector<int,3> vector(ToStart()<<1<<2.5<<3);//warning about conversion double->int

cptG

[31/10/2010 at 23:41:54]

Here's my try, I'm taking advantage of the fact that a function template determines the type of the template parameter from its argument. The only drawback: the numerical type of the array is determined by the first number passed. so arg()=1,2.0,3.0 will fail since it expects integers.
Another caveat is the fact that the initialization-argument needs to be enclosed in brackets, but this is as close as i could get from the top of my head.

Cheers,
Thomas

#include <iostream>
#include <cassert>

namespace detail
{

template <class T,unsigned N>
class vArgT{
public:
    vArgT(vArgT<T,N-1>* left,T arg)
        :left(left)
        ,arg(arg)
    {}

    vArgT<T,N+1> operator,(T arg){
        return vArgT<T,N+1>(this,arg);
    }

    T operator[](unsigned i){
        return i==N?arg:left->operator [](i);
    }

private:
    vArgT<T,N-1>* left;
    T arg;
};

template<class T>
class vArgT<T,0>{
public:
    vArgT<T,0>(T arg)
        :arg(arg)
    {}

    vArgT<T,1> operator,(T arg){
        return vArgT<T,1>(this,arg);
    }

    T operator[](unsigned i){
        assert(i==0);
        return arg;
    }

private:
    T arg;
};

}//detail

class vArg{
public:
    template<class T>
    detail::vArgT<T,0> operator=(T arg){
        return detail::vArgT<T,0>(arg);
    }
};

template <typename T,unsigned N>
class Vector
{
public:
    Vector(detail::vArgT<T,N-1> args)
    {
        for(unsigned i=0;i<N;++i)
            elements[i]=args[i];
    }

    T& operator[](unsigned i){
        return elements[i];
    }

private:
    T elements[N];
};


int main(int argc, char *argv[])
{
    detail::vArgT<double,4> test=(vArg()=3.0,4.0,5.0,6.0,7.0);
    for(int i=0;i<5;++i)
        std::cout<<test[i];

    Vector<double,4> vec((vArg()=1.0,2.0,3.0,4.0));
    for(int i=0;i<4;++i)
        std::cout<<vec[i];
}
(optional)
(optional)
(required, hint)

Links can be added like [this one -> http://www.mr-edd.co.uk], to my homepage.
Phrases and blocks of code can be enclosed in {{{triple braces}}}.
Any HTML markup will be escaped.