// generic column-oriented N-dimensional array of T
struct NDArray(T,int N) {
  T *data;
  uint[N] length;

  // construct array from data and dimensions x0, x1, ...
  static NDArray opCall(T* data, uint x0,...) {
    NDArray res;
    uint[] len = (&x0)[0..N]; // not very portable, use stdarg
    res.length[] = len[];
    if (data is null) {
      uint total = 1;
      for (int k = 0; k < N; ++k) {
	total *= len[k];
      }
      data = new T[total];
    }
    res.data = data;
    return res;
  }

  // private struct for indexing expressions specialized for last dimension
  private struct IndexExpr(int M:1) {
    T* ind_data;
    uint sum;
    uint prod;

    static .NDArray!(T,N).IndexExpr!(1) make(T* data, 
					     uint[] len, 
					     uint sum, 
					     uint prod) {
      .NDArray!(T,N).IndexExpr!(1) res;
      res.ind_data = data;
      res.sum = sum;
      res.prod = prod;
      return res;
    }

    T opIndex(uint n) {
      return ind_data[sum + prod*n];
    }

    void opIndex(uint n,T val) {
      ind_data[sum + prod*n] = val;
    }
  }

  // private struct for indexing expressions for reducing dimension by 1
  // with luck this can be inlined
  private struct IndexExpr(uint M) {
    uint[] ind_length;
    T* ind_data;
    uint sum;
    uint prod;

    static .NDArray!(T,N).IndexExpr!(M) make(T* data, 
					     uint[] len, 
					     uint sum, 
					     uint prod) {
      .NDArray!(T,N).IndexExpr!(M) res;
      res.ind_data = data;
      res.ind_length = len;
      res.sum = sum;
      res.prod = prod;
      return res;
    }

    .NDArray!(T,N).IndexExpr!(M-1) opIndex(uint n) {
      return .NDArray!(T,N).IndexExpr!(M-1).make(ind_data,
						 ind_length[1..M], 
						 prod*n + sum, 
						 prod*ind_length[0]);
    }
  }

  .NDArray!(T,N).IndexExpr!(N-1) opIndex(int n) {
    return .NDArray!(T,N).IndexExpr!(N-1).make(data,length[1..N],n,length[0]);
  }

}


// Two dimensional array is specialized for easier inlining
struct NDArray(T,int N:2) {
  T *data;
  uint[2] length;

  // construct a new matrix
  static NDArray opCall(uint x0, uint x1) {
    return opCall(new T[x0*x1],x0,x1);
  }
  
  static NDArray opCall(T* data, uint x0, uint x1) {
    NDArray res;
    res.length[0] = x0;
    res.length[1] = x1;
    res.data = data;
    return res;
  }

  // private intermediate structure for indexing expressions
  private struct IndexExpr {
    NDArray *x;
    uint m;

    T opIndex(uint n) {
      // bounds-check on length
      return x.data[m + x.length[0]*n];
    }

    void opIndex(uint n,T val) {
      // bounds-check on length
      x.data[m + x.length[0]*n] = val;
    }
  }

  IndexExpr opIndex(int m) {
    IndexExpr res;
    res.x = cast(NDArray*)&data; // no "this" for struct, hmmm
    res.m = m;
    return res;
  }

  // reshape matrix, just for fun
  NDArray reshape(uint x0, uint x1) {
    NDArray res;
    res.length[0] = x0;
    res.length[1] = x1;
    res.data = data;
    return res;
  }

}

template Matrix(T) {
  alias NDArray!(T,2) Matrix;
}