[WIP] te::array<T,int...> discussion

[Remi123]

Preface

This is a blog entry on how to make a variadic multidimensional array class. It will use Template Meta-Programming ( TMP ) and expression templates in the most digestible "slide code" I could.

Introduction to the problem

std::array<typename T, int N> is interesting from the point of view of the C++11 standard library. The Iterator concept was written to mimic pointer, which arrays shares the semantic for iteration, but it wasn't included into the 98 standard.
Maybe it's because it was too "trivial", but we still doesn't have an multidimensional array container, aka array<typename T,int ...Ns>.

Turns out multidimensional arrays are not trivial to implement, as we will see shortly.

An array<T,std::size_t N> example

As a simple reference, let's see a very barebone std::array. There is other functions but those are the main ones.

template<typename T, int N> 
struct array{  // No Ctor because we want array to be an aggregate type.
  using value_type = T;  using size_type = std::size_t; using iterator = T*; using reference = value_type&;
  // Other typedef, but these are the important.
  static constexpr std::size_t size() {return N;}
  value_type buffer[size()]; // The buffer ( private )
  iterator begin() {return &buffer[0];} // Const version cbegin() is almost the same thing
  iterator end() {return &buffer[size()];} // Idem for cend()
  reference operator[](std::size_t index){return &buffer[index];} // Member access operator
  reference at(std::size_t n);// Same as operator[] but we throw if out of range
}; // 

This is nice but let's see what happen when you tries to use variadic to make it multidimensional.

template<typename T, int ... Ns>
struct array{
  // Same typedef as before.
  // static constexpr std::size_t size() {return ...;} // Oups, we need to multiply all Ns...
  // value_type buffer[size()]; // Can't define buffer if size is not working. 
  // iterator begin() {return &buffer[0];} // 
  // iterator end() {return &buffer[size()];} // Idem for cend()
  // reference operator[](std::size_t index); //TMP WARNING
  // reference at(std::size_t n); // TMP WARNING
};

In other words : We need to calculate the size of the array. because anything else depend on it. And the member access operator is a little bit more complicated ( that's a lie, it a lot more complicated ).

Implementing array<T,std::size_t ... Is>::size() ; TMP to the rescue.

Ok so size() can be implemented differently depending on which standard you have on your compiler.

template<std::size_t ...>struct multiply_all // Outside of array struct.
{   static constexpr std::size_t value = 1;};
template<std::size_t I, std::size_t ... Is> struct multiply_all<I,Is...> 
{   static constexpr std::size_t value = I * multiply_all<Is...>::value; };
// In array<T,int...>
static constexpr std::size_t size()const noexcept {return multiply_all<Ns...>::value;} // C++11
static constexpr std::size_t size()const noexcept {return (Ns * ...);}// C++17

Of course, since I have a TMP library, I can simply write this in C++11 with using namespace te;

static constexpr std::size_t size(){return eval_pipe_<input_<i<1>>,multiply_<I<Ns>>...>::value;} // C++11

The rest of the function except member access operator[] can be simply written like this :

template<typename T, int ... Ns>
struct array{
  // Same typedef as before.
  static constexpr std::size_t size(){return eval_pipe_<input_<i<1>>,multiply_<I<Ns>>...>::value;} 
  value_type buffer[size()];  
  iterator begin() {return &buffer[0];} 
  iterator end() {return &buffer[size()];} 
  // reference operator[](std::size_t index); // TMP WARNING
  // template<typename ... Size_t>reference at(Size_t ... ns); // TMP WARNING
};

Calculating the strides

As of now, we have a working buffer that can be iterated over, using the standard library functions. But that's boring and I know you are there for the member access operator. The problem is that variadic expansion doesn't work with the square bracket like this :

reference operator[](const std::size_t ... indices) {return buffer[indices]...;}

It doesn't expand into buffer[N0][N1][N2]... or even buffer[N0,N1,N2,...], it just crash. There exist a proposals to "make this work" by also deprecating the use of overloaded comma iterator in square bracket operator. (Overloading the comma operator is always a bad idea). The proposed mdspan and other library use the call operator with variadic parameter to simulate this, and boost.MultiArray use a temporary class that iterate between each [index]. The boost solution is the correct one, but share more similarities with Expression Template than normal code.

Let's program the at(std::size_t ... n) variadic templated function. But first we need some helper constexpr function similar to size(), or rather we need create arrays with information at compile time.

template<typename ... Ts>
    struct array_value_constant
    {
        constexpr std::size_t operator[](const std::size_t& n)const{
            constexpr std::size_t val[] = {(std::size_t)Ts::value...};
            return val[n];}
    };
template<typename T,std::size_t...Ns>
struct array{ //..
  static constexpr inline std::size_t dimension_size()noexcept   {return sizeof...(Is);}
  static constexpr std::size_t shapes_c[] = {Is...};
  static constexpr inline const std::size_t shapes(const std::size_t n)noexcept {   return shapes_c[n];}
  static constexpr inline const std::size_t strides(const std::size_t n)noexcept 
  { constexpr auto array_value = te::eval_pipe_<
      te::input_<std::integral_constant<std::size_t,size()>,std::integral_constant<std::size_t,Is>...>
      ,te::fold_left_list_<te::divide_<>>
      ,te::wrap_<array_value_constant>>;
    return array_value[n];  }
  //...  };

1. dimension_size() is just returning the number of dimension that the array possess. Nothing too complicated.
2. The constexpr function const std::size_t* shapes() (notes: name come from boost.multiarray) returns an array of size sizeof...(Is) that is the same 'shapes' as the variadic to help iterating over them. (e.g. array<int,4,5,6>::shapes() == {4,5,6}). We use a constexpr variable to form it, but this technique have limit if we want more functionality for the next function.
3. This one is more complicated. We have to create an array of the strides of the array, basically the numbers of that allow use to calculate multiple dimension index into a single dimension index. For example, array<int,2,5,4>{...}[1][2][3] has a size of 40 (2 * 5 * 4), and to calculate an index in a single dimension, we want to do buffer[1*20+2*4+3*1] == buffer[31]. The strides are the number {20,4,1} for this array. Basically, we want to take the size of the array, then recursively divide it by each shapes and save each intermediate result into an array. 40/2/5/4 ->{40/2==20, 20/5==4, 4/4==1} -> {20,4,1}. It should always end with a value of 1. This is technically a fold expression where you append the result to a vector, but I'm in C++14 so no luck I need to do it myself.
4. To create an C array from that, I use the metafunction array_value_constant and ship it all the integral_constant required (this is the wrap_<array_value_constant> part). Those integral_constant come from input_<i<size()>,i<Is>...> and finally we pass the meta-function te::fold_left_list_<te::divide_<>> which does all the hard work in a generic way (for my library at least). This result in an array_value_constant<Strides_c...>. To bypass some linker issues, the member access operator create a constexpr array, and return the value at the correct index.

That's a lot for simply creating a C array but it work correctly.

Fortunately the at(std::size_t ... n) is rather simple to implement once we have this.

    template<typename ... Size_t>
    constexpr reference at(Size_t ... index )noexcept{
        static_assert(sizeof...(Size_t) == dimension_size(),"Number of argument must equal the number of dimensions");
        std::size_t indices[] = {(std::size_t)index...};
        std::size_t result = 0;
        for(std::size_t i = 0; i < dimension_size();++i )
            if(indices[i] < shapes(i) {result += indices[i] * strides(i);} 
            else {throw std::out_of_range("array member access is higher than the maximum of its dimension");}
        return buffer[result];
    }

This isn't so bad in my opinion. We throw std::out_of_range if one of the index is higher than its respective dimension, otherwise we increase the result by the product of the index and stride at the same index.

Implementing typical member access operator[] ( my_array_3dim[1][2][3];)

Next, the clou du spectacle, the member access operator[]. For this, we need to go into "expression template", which is a fancy way of creating new behavior by creating "temporary class" in-between expressions. In the case of the expression my_array[1][0][2] from an array<int,2,3,4>, each [i] except the last one returns an instance of a class that keep a reference to the buffer, and modify an resulting index depending on its strides. The last one returns the reference to the buffer data at the correct index.

template<typename T, int ... Is>
struct array{
// ...
    // Operator []
    template <bool Pred = sizeof...(Is) == 1>
    constexpr auto operator[](std::size_t n) const noexcept
        -> std::enable_if_t<Pred, const_reference>
    {        return buffer[n];    }
private:
    // Member access private class to make it work with arr[0][1][2]
    template<std::size_t I,std::size_t N = array::dimension_size()-1>
    struct expr_tmp_mem_acc_op{
        const value_type* data;
        std::size_t index;
        constexpr expr_tmp_mem_acc_op<I+1,N> operator[](std::size_t n)noexcept
        {return {data,index+n*strides(I)};}
        constexpr const expr_tmp_mem_acc_op<I+1,N> operator[](std::size_t n)const noexcept
        {return {data,index+n*strides(I)};}
    };
    template<std::size_t N>
    struct expr_tmp_mem_acc_op<N,N>
    {   const value_type* data;
        std::size_t index;
        constexpr reference operator[](std::size_t n)noexcept
        {return const_cast<reference>(data[index + n*strides(N)]);}
        constexpr const_reference operator[](std::size_t n)const noexcept
        {return data[index + n*strides(N)];}
    };
public:
    template <bool Pred = (sizeof...(Is) > 1)>
    constexpr auto operator[](std::size_t n) noexcept
        -> std::enable_if_t<Pred, expr_tmp_mem_acc_op<1>>
    { return {buffer,n * strides(0)};}
    template <bool Pred = (sizeof...(Is) > 1)>
    constexpr auto operator[](std::size_t n) const noexcept
        -> std::enable_if_t<Pred, const expr_tmp_mem_acc_op<1>>
    { return {buffer,n * strides(0)};}
// ..
};

Ok this is a lot but it's the clearer I could do. The first operator[] is only available through SFINAE if the number of dimension is one and it does the typical array access. But for the case we have more than one dimension, I need the template class expr_tmp_mem_acc_op<std::size_t, std::size_t = array::dimension_size()> to "count" the number of square bracket until the second to last, where in this case the class expr_tmp_mem_acc_op<N,N> is specialized to return the actual reference. Each square bracket updates the calculation of the index, and the last one return the reference to the value.

Some of you may be concerned about the const_cast in {return const_cast<reference>(data[index + n*strides(N)]);}, with good reason as every const_cast must be justified. Fortunately this path is only taken if and only if the instance of the array is const, as you can see with const expr_tmp_mem_acc_op<1> operator[](std::size_t)const noexcept when const and expr_tmp_mem_acc_op<1> operator[](std::size_t) noexcept when non-const. My reasoning is that every expr_tmp_mem_acc store the pointer to the buffer as a const pointer so instead of propagating constness, I have const as default and const_cast if non-const.

Empty case

The only thing left to write about is the strange case of array<T>, array<T,0> and array<T,0,...,0>. After much consideration, I will not letarray<T> compile. as it's an array with no dimension and I can't reason about. The other two have dimensions, but with the size of zero. You can write constexpr auto a = array<int,0,0>{};, but the constructor can't take arguments. They will have an unitialized memory.

Conclusion

There is a couple of things I did not show in this blog, like most of the function have a const equivalent, the strides meta-functions is slightly more complex to avoid division by zero, but overall what you see here is mostly the same as what is in my library. I hope you found this blog interesting, since it touches on the subject of meta-programming, expression templates and simple memory access.