Variadic expansion to access multi-dimensional std::array

Question

Suppose the following multi_array structure:

template <typename type, std::size_t... sizes>
struct multi_array
{
  using storage_type    = typename storage_type<type, sizes...>::type;
  using value_type      = type;
  using size_type       = std::array<std::size_t   , sizeof...(sizes)>;
  using difference_type = std::array<std::ptrdiff_t, sizeof...(sizes)>;
  using reference       = value_type&;
  using const_reference = const value_type&;
  // ...

  storage_type data_  {};
  size_type    sizes_ {sizes...};
}

// Example usage: 
// multi_array<float, 3, 4, 5> three_dimensional_float_array;
// storage_type will be std::array<std::array<std::array<float, 5>, 4>, 3>
// size_type    will be std::array<std::size_t, 3>

where:

// Helpers to create a std::array<std::array<...>> from a pack expansion.
template <typename  type,                   std::size_t... sizes>
struct storage_type;
template <typename _type, std::size_t size, std::size_t... sizes>
struct storage_type<_type, size, sizes...>
{
  using type = std::array<typename storage_type<_type, sizes...>::type, size>;
};
template <typename _type>
struct storage_type<_type>
{
  using type = _type;
};

I am now trying to implement the .at() function:

  [[nodiscard]]
  constexpr reference       at   (const size_type& position)
  {
    // This has to be:
    // data_[position[0]][position[1]][position[2]]...;
  }

Fold expressions do not work in this case so a recursive solution is necessary.

I think this is the way, but can't seem to come up with the right answer:

[[nodiscard]]
constexpr reference at(const size_type& position)
{
  return access_type<value_type, sizeof...(sizes)>::at(data_, position);
}

template <typename type, std::size_t size, std::size_t index = 0>
struct access_type
{
  template <typename array_type>
  auto at (const array_type& array, const std::array<std::size_t, size>& position, type& current) // Problem: Reference is not a type& for intermediates.
  {
    if constexpr (index == 0)
      current = &array;

    if constexpr (index + 1 != size)
    {
      return access_type::at<type, size, index + 1>(array, position, current);
    }
  }
};

Answer 1

Something along these lines, perhaps:

template <size_t level, size_t max_level>
struct Access {
  template<typename Store>
  static reference At(Store& store, const size_type& position) {
    return Access<level + 1, max_level>::At(
        store[position[level]], position);
  }
};

template <size_t level>
struct Access<level, level> {
  template<typename Store>
  static reference At(Store& store, const size_type&) {
    return store;
  }
};


reference at(const size_type& position)
{
  return Access<0, sizeof...(sizes)>::At(data_, position);
}

Answer 2

Here is my take:

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);
    }
}

It transforms tuple (x,y,z,...) into 1-D index x*N + y*K + ..., where N, K, ... are chosen accordingly to the total container dimensions.

And used in my own multi-dimension cache-friendly container:

    /**
     * @brief Returns a reference to the element at coordinates.
     * @param coordinates Coordinates of the element to return
     *
     * No bounds checking is performed; if @c coordinates are outside od
     * the matrix dimensions, the behavior is undefined.
     */
    template<class... Args>
    T const& operator()(Args... coordinates) const
    { return _data[_details::coordinates_to_index(dimensions, std::array{coordinates...})]; }

Source : yscialom/matrix/matrix.hpp on my GitHub (PR welcomed)