mdspan.hpp

template <typename T, size_t N>
struct mdspan {
  T *data;
  std::array<size_t, N> extents;

  template <typename... OtherIndexTypes>
    requires(std::is_convertible_v<OtherIndexTypes, size_t> && ...) &&
                (sizeof...(OtherIndexTypes) == N)
  constexpr explicit mdspan(T *data, OtherIndexTypes... _extents)
      : data(data), extents{} {
    auto idx = size_t{0};
    (
        [&](OtherIndexTypes extent) {
          extents[idx++] = static_cast<size_t>(extent);
        }(_extents),
        ...);
  }

  template <typename... OtherIndexTypes>
    requires(std::is_convertible_v<OtherIndexTypes, size_t> && ...) &&
            (sizeof...(OtherIndexTypes) == N)
  constexpr auto operator()(OtherIndexTypes... indices) -> T & {
    auto real_idx = size_t{0};
    auto ext_idx = size_t{0};

    (
        [&](OtherIndexTypes idx) {
          real_idx *= extents[ext_idx++];
          real_idx += static_cast<size_t>(idx);
        }(indices),
        ...);

    return data[real_idx];
  }

#if __cpp_multidimensional_subscript >= 202110L
  template <typename... OtherIndexTypes>
    requires(std::is_convertible_v<OtherIndexTypes, size_t> && ...) &&
            (sizeof...(OtherIndexTypes) == N)
  constexpr auto operator[](OtherIndexTypes... indices) -> T & {
    return operator()(indices...);
  }
#endif
};

template <typename T, typename... Integrals>
  requires(sizeof...(Integrals) > 0)
explicit mdspan(T *, Integrals...) -> mdspan<T, sizeof...(Integrals)>;