#include <cassert>
#include <iostream>
#include <vector>

template <typename ELEM>
class NDArrayT {
  private:
    // dimensions
    std::vector<size_t> _dims;
    // data
    std::vector<ELEM> _data;
    
  public:
    NDArrayT(const std::vector<size_t> &dims):
      _dims(dims)
    {
      size_t size = _dims.empty() ? 0 : 1;
      for (size_t dim : _dims) size *= dim;
      _data.resize(size);
    }
    NDArrayT(
      const std::vector<size_t> &dims,
      const std::vector<ELEM> &data):
      NDArrayT<ELEM>(dims)
    {
      assert(_data.size() == data.size());
      std::copy(data.begin(), data.end(), _data.begin());
    }

    ELEM& operator[](const std::vector<size_t> &indices)
    {
      size_t i = 0, j = 0;
      for (size_t n = _dims.size(); j < n; ++j) {
        i *= _dims[j]; i += indices[j];
      }
      return _data[i];
    }

    const ELEM& operator[](const std::vector<size_t> &indices) const
    {
      size_t i = 0, j = 0;
      for (size_t n = _dims.size(); j < n; ++j) {
        i *= _dims[j]; i += indices[j];
      }
      return _data[i];
    }
};

using namespace std;

ostream& operator<<(ostream &out, const vector<size_t> &values)
{
  const char *sep = "";
  for (size_t value : values) {
    out << sep << value; sep = ", ";
  }
  return out;
}

bool inc(vector<size_t> &indices, const vector<size_t> &dims)
{
  for (size_t i = indices.size(); i--;) {
    if (++indices[i] < dims[i]) return false;
    indices[i] = 0;
  }
  return true; // overflow
}

int main()
{
  // build sample data
  vector<double> data(2 * 3 * 4);
  for (size_t i = data.size(); i--;) data[i] = (double)i;
  // build sample array
  typedef NDArrayT<double> NDArrayDouble;
  const vector<size_t> dims = { 2, 3, 4 };
  NDArrayDouble a(dims, data);
  // print sample array (check subscript)
  vector<size_t> indices(dims.size(), 0);
  do {
    cout << "a[" << indices << "]: " << a[indices] << endl;
  } while (!inc(indices, dims));
  // done
  return 0;
}