# A basic multi-dimensional array

I have been lurking on the C++ side of stack-overflow only long enough to know there are a lot of beginners and intermediate programmers baffled by multi-dimensional arrays. I've seen a lot of monstrosities : three-star programmers, cache-unfriendly implementations, vectors of vectors of vectors ..., etc.

I thought to provide them with a basic multi-dimensional array:

• open for extensions,
• fit for production and experimentation,
• working without allocations,
• and cache-friendly.

Before I do so, I'd like some advice from the community.

1. What can I improve in the provided features of this class?
2. What new feature could one add, and why would it be useful?
#include <array>
#include <numeric>

namespace ysc
{
namespace _details
{
template<class InputIt, class OutputIt>
OutputIt partial_product(InputIt first, InputIt last, OutputIt output)
{ *output++ = 1; return partial_sum(first, last, output, std::multiplies<>{}); }

// cache-friendly:
// neighbor objects within the right-most coordinate are neighbors in memory
template<class TDim, class TCoord>
auto coordinates_to_index(TDim const& dimensions, TCoord const& coords)
{
std::array<std::size_t, dimensions.size()> dimension_product;
using std::crbegin, std::crend, std::prev;
partial_product(crbegin(dimensions), prev(crend(dimensions)), begin(dimension_product));
return std::inner_product(cbegin(dimension_product), cend(dimension_product), crbegin(coords), 0);
}
}

constexpr struct matrix_zero_t {} matrix_zero;

template<class T, std::size_t... Dimensions>
class matrix
{
template<class, std::size_t...> friend class matrix;

public:
static constexpr std::size_t order      = sizeof...(Dimensions);
static constexpr std::array  dimensions = { Dimensions... };

private:
std::array<T, (Dimensions * ...)> _data;

public:
friend void swap(matrix& lhs, matrix& rhs)
{
using std::swap;
swap(lhs._data, rhs._data);
}

public:
matrix()                          = default;
matrix(matrix&& other)            = default;
matrix& operator=(matrix&& other) = default;

matrix(matrix_zero_t) : _data({}) {}

template<class U>
matrix(matrix<U, Dimensions...> const& other) { std::copy(cbegin(other._data), cend(other._data), begin(_data)); }

template<class U>
matrix& operator=(matrix<U, Dimensions...> const& other)
{
matrix o{other};
swap(*this, o);
return *this;
}

public:
template<class... Args>
T const& operator()(Args... coordinates) const
{ return _data[_details::coordinates_to_index(dimensions, std::array{coordinates...})]; }

template<class... Args>
T& operator()(Args... coordinates)
{ return _data[_details::coordinates_to_index(dimensions, std::array{coordinates...})]; }
};
} // namespace ysc


This is a C++17 implementation; this itself is not set in stone.

• How would this be used? – Mast Oct 26 '18 at 8:45
• @Mast I've included a small example in the usage demo I link under the code. This would be used as a class with value semantics: matrix<...> m = matrix_zero; m(x,y,z) = some_value; m(a,b,c) = m(x,y,z) + m(y,z,x);. – YSC Oct 26 '18 at 8:50
• @YSC error: cannot bind non-const lvalue reference of type 'ysc::matrix<float, 2>&' to an rvalue of type 'ysc::matrix<float, 2>'. Apparently not. And non-converting copy construction / assignment don't compile either (for other reasons). – user673679 Oct 26 '18 at 9:58
• @user673679 You're completely right. Shame on me! (also: fixed) – YSC Oct 26 '18 at 10:03
• While I like the idea behind this class, I don't think it's better than some nested std::arrays (for example like this rough sketch) with its current feature set. // In my experience, when working with multidimensional arrays you often want some way to refer to some "subslice" of the covered space. That would be a nice feature ;) – hoffmale Oct 26 '18 at 20:46

Here are some suggestions:

1. Your coordinates_to_index should better return a std::size_t instead of int because int may not be able to hold the required values.

return std::inner_product(cbegin(dimension_product), cend(dimension_product), crbegin(coords),
std::size_t(0)); // note: not int

2. You don't need to compute the strides and store them. You can use Horner's rule directly.

std::size_t result = 0;
for (std::size_t i = 0; i < N; ++i)
result = result * dimensions[i] + coords[i];
return result;

3. Your calculation of total size fails to consider the case where order == 0. Instead of

(Dimensions * ...)


Use

(Dimensions * ... * static_cast<std::size_t>(1))


Moreover, I would suggest exposing the total size, along with the order and dimensions. Like this:

static constexpr std::size_t size = (Dimensions * ... * std::size_t(1));


And use it in place of ad-hoc computations further in code.

4. Your matrix is not Copyable. You declare move operations but not copy operations (the template doesn't make a difference), and the compiler synthesizes deleted copy operations which take precedence over the copying template. You should either default the copy operations explicitly (and maybe destructor as well), or omit the redundant default declarations directly.

5. Why use :_data({}) when :_data{} is more readable and less verbose?

6. Consider passing by value instead of making a copy yourself. This enables moving. And this is likely to be optimized better.

template <class U>
matrix& operator=(matrix<U, Dimensions...> other) // pass by value
{
swap(*this, other);
return *this;
}

7. Consider adding out-of-range detection for operator() in some form. It should be an opt-out feature, not opt-in.