//
// rarray - Runtime arrays: template classes for pointer-based,
//          runtime dimensionalized, multi-dimensional
//          arrays.  Documentation in rarraydoc.pdf
//
// Copyright (c) 2013-2020  Ramses van Zon
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
// THE SOFTWARE.
//

#ifndef RARRAY_H
#define RARRAY_H

#if __cplusplus >= 201103L

#include <stdexcept>
#include <cstdlib>
#include <utility>
#include <string>
#include <cmath>

#ifdef RA_BOUNDSCHECK
#include <string>
#define RA_CHECKORSAY(a, b) if (not(a)) throw std::out_of_range(std::string(b) + " in function " + std::string(__PRETTY_FUNCTION__))
#else
#define RA_CHECKORSAY(a, b) 
#endif

#if not defined(RA_INLINE)
# if defined(__INTEL_COMPILER)
#   define RA_INLINE  __forceinline
# elif defined(__GNUC__)
#   define RA_INLINE __attribute__((always_inline)) inline
# else
#   define RA_INLINE inline
# endif
#endif

#define RA_INLINEF RA_INLINE
#define RA_INLINE_ inline



#ifndef SHAREDBUFFERH
#define SHAREDBUFFERH

#include <cstddef>
#include <algorithm>
#include <stdexcept>
#include <memory>
#include <cstring>
#include <atomic>

int test_shared_buffer_main(); // for testing

/***************************************************************************/

namespace ra {
static char _buffername = 'A'; 

/***************************************************************************/

template<class T>
class shared_buffer
{
  public:

    typedef ssize_t size_type;

    shared_buffer();
    explicit shared_buffer(size_type size);
    shared_buffer(size_type size, T* data);
    shared_buffer(const shared_buffer& other);
    shared_buffer(shared_buffer&& from);
    
    shared_buffer& operator=(const shared_buffer& other);  // copy
    void operator=(shared_buffer&& from);                  // move

    ~shared_buffer();

    const T& operator[](size_type index) const;
    T& operator[](size_type index);
    
    const T& at(size_type index) const;
    T& at(size_type index);

    shared_buffer<T> slice(size_type from, size_type to);

    size_type size() const;

    shared_buffer<T> copy() const;

    typedef T* iterator;
    typedef const T* const_iterator;
    iterator begin();
    iterator end();
    const_iterator begin() const;
    const_iterator end() const;
    const_iterator cbegin() const;
    const_iterator cend() const;

    typedef std::reverse_iterator<iterator> reverse_iterator;
    typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
    reverse_iterator rbegin();
    reverse_iterator rend();
    const_reverse_iterator rbegin() const;
    const_reverse_iterator rend() const;
    const_reverse_iterator crbegin() const;
    const_reverse_iterator crend() const;

    void resize(size_type newsize, bool keep_content = false);
    
    void assign(const T& value);
    void assign(std::initializer_list<T> ilist);
    template<class InputIt>
    void assign_iter(InputIt first, InputIt last);
   
  private:
    
    T*        data_;
    T*        orig_;
    size_type size_;
    std::atomic<int>* refs_;

    friend int ::test_shared_buffer_main(); // for testing

    void uninit();
    void incref();
    void decref();
};

/***************************************************************************/
#define RA_ALIGNMENT_IN_BYTES 64

template<class T>
void malign(void*& orig, void*& place, size_t size, size_t bytealignment)
{
    size_t paddedsizeinbytes = size*sizeof(T) + bytealignment - 1;
    void* some_orig = malloc(paddedsizeinbytes);
    if (some_orig) {        
        void* some_place = some_orig;
        if (std::align(bytealignment,
                       size*sizeof(T),
                       some_place,
                       paddedsizeinbytes)) {
            orig = some_orig;
            place = some_place;
        } else {
            orig = nullptr;
            place = nullptr;
        }
    }
}

template<class T>
shared_buffer<T>::shared_buffer()
{
    uninit();
}

template<class T>
shared_buffer<T>::shared_buffer(size_type size)
    : data_(new T[size]), orig_(data_), size_(size), refs_(new std::atomic<int>(1))
{
}

template<class T>
shared_buffer<T>::shared_buffer(size_type size, T* data)
  : data_(data), orig_(nullptr), size_(size), refs_(nullptr)
{
}

template<class T>
shared_buffer<T>::shared_buffer(const shared_buffer<T>& other)
  : data_(other.data_), orig_(other.orig_), size_(other.size_), refs_(other.refs_)
{
    incref();
}

template<class T>
shared_buffer<T>::shared_buffer(shared_buffer<T>&& from)
  : data_(from.data_), orig_(from.orig_), size_(from.size_), refs_(from.refs_)
{
    from.uninit();
}

template<class T>
shared_buffer<T>& shared_buffer<T>::operator=(const shared_buffer<T>& other)
{
    if (this != &other) {
        decref();
        data_ = other.data_;
        orig_ = other.orig_;
        size_ = other.size_;
        refs_ = other.refs_;
        incref();
    }
    return *this;
}

template<class T>
void shared_buffer<T>::operator=(shared_buffer<T>&& from)
{
    decref();
    data_ = from.data_;
    orig_ = from.orig_;
    size_ = from.size_;
    refs_ = from.refs_;
    from.uninit();
}

template<class T>
const T& shared_buffer<T>::operator[](size_type index) const
{
    return data_[index];
}

template<class T>
T& shared_buffer<T>::operator[](size_type index)
{
    return data_[index];
}

template<class T>
const T& shared_buffer<T>::at(size_type index) const
{
    if (index < 0 or index > size_)
        throw std::out_of_range("shared_buffer::at");
    return data_[index];
}

template<class T>
T& shared_buffer<T>::at(size_type index)
{
    if (index < 0 or index > size_)
        throw std::out_of_range("shared_buffer::at");
    return data_[index];
}

template<class T>
shared_buffer<T> shared_buffer<T>::slice(size_type from, size_type to)
{
    if (from < 0 or to < 0 or from > size_ or to > size_)
        throw std::out_of_range("shared_buffer::slice");
    shared_buffer<T> result(*this);
    result.data_ += from;
    if (from <= to)
        result.size_ = to - from;
    else
        result.size_ = 0;
    return result;
}

template<class T>
shared_buffer<T>::~shared_buffer()
{
    decref();
}

template<class T>
typename shared_buffer<T>::size_type shared_buffer<T>::size() const
{
    return size_;
}

template<class T>
shared_buffer<T> shared_buffer<T>::copy() const    
{
    shared_buffer<T> result(size_);
    std::copy(cbegin(), cend(), result.begin());
    return result;
}


template<class T>
typename shared_buffer<T>::iterator shared_buffer<T>::begin()
{
    return data_;
}

template<class T>
typename shared_buffer<T>::iterator shared_buffer<T>::end()    
{
    return data_+size_;
}

template<class T>
typename shared_buffer<T>::const_iterator shared_buffer<T>::cbegin() const
{
    return data_;
}

template<class T>
typename shared_buffer<T>::const_iterator shared_buffer<T>::cend() const
{
    return data_+size_;
}

template<class T>
typename shared_buffer<T>::const_iterator shared_buffer<T>::begin() const
{
    return data_;
}

template<class T>
typename shared_buffer<T>::const_iterator shared_buffer<T>::end() const
{
    return data_+size_;
}


template<class T>
typename shared_buffer<T>::reverse_iterator shared_buffer<T>::rbegin()
{
    return std::reverse_iterator<iterator>(data_+size_);
}

template<class T>
typename shared_buffer<T>::reverse_iterator shared_buffer<T>::rend()
{
    return std::reverse_iterator<iterator>(data_);
}

template<class T>
typename shared_buffer<T>::const_reverse_iterator shared_buffer<T>::crbegin() const
{
    return std::reverse_iterator<const_iterator>(data_+size_);
}

template<class T>
typename shared_buffer<T>::const_reverse_iterator shared_buffer<T>::crend() const
{
    return std::reverse_iterator<const_iterator>(data_);
}

template<class T>
typename shared_buffer<T>::const_reverse_iterator shared_buffer<T>::rbegin() const
{
    return std::reverse_iterator<const_iterator>(data_+size_);
}

template<class T>
typename shared_buffer<T>::const_reverse_iterator shared_buffer<T>::rend() const
{
    return std::reverse_iterator<const_iterator>(data_);
}

template<class T>
void shared_buffer<T>::resize(size_type newsize, bool keep_content)
{
    if ( (newsize<size_) and (refs_) and (*refs_)==1) {
        size_ = newsize; // That's all? Yep!
    } else {
        T* newdata = new T[newsize];
        T* neworig = newdata;
        std::atomic<int>* newrefs = new std::atomic<int>(1);
        if (keep_content) {
            size_t n = ((size_<newsize)?size_:newsize);
            for (size_t i=0; i<n; i++)
                newdata[i] = data_[i];
        }
        decref();
        data_ = newdata;
        orig_ = neworig;
        size_ = newsize;
        refs_ = newrefs;
    }
}

template<class T>
void shared_buffer<T>::assign(const T& value)
{
    for (int i=0;i<size_;i++)
        data_[i] = value;
}

template<class T>
template<class InputIt>
void shared_buffer<T>::assign_iter(InputIt first, InputIt last)
{
    resize(last-first);
    T* data = data_;
    for (InputIt it = first; it != last; it++)
        *(data++) = *it;
}

template<class T>
void shared_buffer<T>::assign(std::initializer_list<T> ilist)
{
    assign_iter(ilist.begin(), ilist.end());
}

/***************************************************************************/

template<class T>
void shared_buffer<T>::uninit() {
    data_ = nullptr;
    orig_ = nullptr;
    size_ = 0;
    refs_ = nullptr;
}

template<class T>
void shared_buffer<T>::incref() {
    if (refs_)
        (*refs_)++;
}

template<class T>
void shared_buffer<T>::decref() {
    if (refs_) {
        if (--(*refs_) == 0) {
            delete[] orig_;
            delete refs_;
            uninit();
        }
    }
}

/***************************************************************************/
}
#endif


#ifndef SHARED_SHAPEH
#define SHARED_SHAPEH

// 
//#include "offsets.h"
#ifndef OFFSETSH
#define OFFSETSH

#include <cstdlib>
#include <vector>
#include <type_traits>

/***************************************************************************/

class Offsets
{
  public:

    Offsets(const std::vector<ssize_t>& dims);

    template<class T>  void*** apply_offsets(T* data) const;

    ssize_t get_num_data_offsets() const;
    ssize_t get_num_offsets() const;
    ssize_t get_rank() const;
    
  private:

    ssize_t rank_;
    std::vector<unsigned long long> offsets_;
    ssize_t ndataoffsets_;

    friend int test_offsets_main(); 
};

/***************************************************************************/

inline Offsets::Offsets(const std::vector<ssize_t>& dims)
{
    std::vector<ssize_t> extent(dims);
    rank_ = extent.size();
    ssize_t noffsets = 0;
    ndataoffsets_ = 0;
    if (rank_ > 0) {
        ndataoffsets_ = 1;
        for (ssize_t i = rank_ - 1; i--; )
            noffsets = extent[i]*(1 + noffsets);
        for (ssize_t i = 0 ; i < rank_-1; i++ ) 
            ndataoffsets_ *= extent[i];
        offsets_.reserve(noffsets);
        offsets_.resize(noffsets);
        if (noffsets > 1) { // catches rank_ == 1 and zero-dimensioned arrays
            ssize_t offsetnum = 0;
            ssize_t extenttot = extent[0];
            for (ssize_t i = 1; i < rank_ - 1; i++) {
                for (ssize_t j = 0; j < extenttot; j++) 
                    offsets_[offsetnum+j] = offsetnum + extenttot + j*extent[i];
                offsetnum += extenttot;
                extenttot *= extent[i];
            }
            ndataoffsets_ = extenttot;
            for (ssize_t j = 0; j < ndataoffsets_; j++) 
               offsets_[offsetnum + j] = j*extent[rank_ - 1];
        }
    }
}

template<class T>
inline void*** Offsets::apply_offsets(T* data) const
{
    ssize_t noffsets = offsets_.size();
    if (ndataoffsets_ == 0 && noffsets == 0)
        return nullptr;
    else if (ndataoffsets_ == 1 && noffsets == 0) // that happens only for rank==1
        return (void***)(data);
    else {
        void*** offsets = new void**[noffsets];
        ssize_t i = 0;
        for (;i < noffsets - ndataoffsets_; i++)            
            offsets[i] = (void**)offsets + offsets_[i];
        for (;i < noffsets; i++)            
            offsets[i] = reinterpret_cast<void**>(
                              const_cast<typename std::remove_const<T>::type*>(data)
                              + offsets_[i]);
        return offsets;
    }
}

inline ssize_t Offsets::get_num_data_offsets() const
{
    return ndataoffsets_;
}

inline ssize_t Offsets::get_num_offsets() const
{
    return offsets_.size();
}

inline ssize_t Offsets::get_rank() const
{
    return rank_;
}

/***************************************************************************/

#endif

//end of #include "offsets.h"
#include <array>
#include <stdexcept>
#include <atomic>

int test_shared_shape_main();

/***************************************************************************/

namespace ra {

/***************************************************************************/


template<typename T,int R> 
struct PointerArray {
    typedef typename PointerArray<T,R-1>::type const*    type;         // const shape, recursive
    typedef typename PointerArray<T,R-1>::noconst_type*  noconst_type; // non-const variant
};
template<typename T> 
struct PointerArray<T,1> { // We end the recursion by specifically defining the R=1 case
    typedef T* type;     // no const as this would express that the elements are constant
    typedef T* noconst_type; 
};
template<typename T> 
struct PointerArray<T,0> { // We end the recursion also by specifically defining the R=0 case
    typedef T& type;     // no const as this would express that the elements are constant
    typedef T& noconst_type; 
};
/***************************************************************************/

template<class T, int R>
class shared_shape
{
  public:

    typedef typename PointerArray<T,R>::type ptrs_type;
    typedef ssize_t size_type;

    shared_shape();                                         // non-functional shape
    shared_shape(const std::array<size_type,R>&extent, T*data);// construct shape 
    shared_shape(const shared_shape& other);                // copy constructor
    shared_shape(shared_shape&& other);                     // move constructor

    shared_shape& operator=(const shared_shape& other);     // copy
    void operator=(shared_shape&& other);                   // move

    ~shared_shape();

    shared_shape copy() const; 

    void relocate(T* newdata);                              // let shape point to other data block
    void reshape(const std::array<size_type,R>&extent);     // change shape (not the data)

    ptrs_type ptrs() const;                                 // get pointer-to-pointer structure
    T* data() const;                                        // get pointer to the data
    size_type size() const;                                 // get total number of elements
    size_type extent(int i) const;                          // get extent in dimension i
    const std::array<size_type,R>& extent() const;          // get total extent array
    
    shared_shape<T,R-1> at(size_type r) const;
    
  private:

    std::array<size_type,R> extent_;
    ptrs_type ptrs_;
    std::atomic<int>* refs_;
    void***   orig_;
    size_type noffsets_;
    size_type ndataoffsets_;

    void uninit();
    void incref() const;
    void decref();
    void copy_before_write();

    template<class U, int S> friend class shared_shape; // for "at"

    friend int ::test_shared_shape_main();
};

template<typename T> 
class shared_shape<T,0> {
  public:
    typedef T ptrs_type;
    typedef ssize_t size_type;
  private:
    std::array<size_type,0> extent_;
    ptrs_type ptrs_;
    std::atomic<int>* refs_;
    void***   orig_;
    size_type noffsets_;
    size_type ndataoffsets_;
    void incref() {}
    template<class U, int S> friend class shared_shape; // for "at"
};
    
/***************************************************************************/

template<class T, int R>
shared_shape<T,R>::shared_shape()
{
    uninit();
}

template<class T, int R>
shared_shape<T,R>::shared_shape(const std::array<size_type,R>&extent, T*data)
  : extent_(extent), ptrs_(nullptr), refs_(nullptr), orig_(nullptr)
{
    Offsets P({extent_.begin(),extent_.end()});
    orig_ = P.apply_offsets(data);
    ptrs_ = reinterpret_cast<ptrs_type>(orig_);
    noffsets_ = P.get_num_offsets();
    ndataoffsets_ = P.get_num_data_offsets();
    if (R>1) {
        refs_ = new std::atomic<int>(1);
    } else {
        orig_ = nullptr;
    }
}

template<class T, int R>
shared_shape<T,R>::shared_shape(const shared_shape& other)
  : extent_(other.extent_),
    ptrs_(other.ptrs_),
    refs_(other.refs_),
    orig_(other.orig_),
    noffsets_(other.noffsets_),
    ndataoffsets_(other.ndataoffsets_)
{
    incref();
}

template<class T, int R>
shared_shape<T,R>::shared_shape(shared_shape&& other)
  : extent_(other.extent_),
    ptrs_(other.ptrs_),
    refs_(other.refs_),
    orig_(other.orig_),
    noffsets_(other.noffsets_),
    ndataoffsets_(other.ndataoffsets_)
{
    other.uninit();
}

template<class T, int R>
shared_shape<T,R>& shared_shape<T,R>::operator=(const shared_shape& other)
{
    if (this != &other) {
        decref();
        extent_       = other.extent_;
        ptrs_         = other.ptrs_;
        refs_         = other.refs_;
        orig_         = other.orig_;
        noffsets_     = other.noffsets_;
        ndataoffsets_ = other.ndataoffsets_;
        incref();
    }
    return *this;
}

template<class T, int R>
void shared_shape<T,R>::operator=(shared_shape&& other)
{
    decref();
    extent_       = other.extent_;
    ptrs_         = other.ptrs_;
    refs_         = other.refs_;
    orig_         = other.orig_;
    noffsets_     = other.noffsets_;
    ndataoffsets_ = other.ndataoffsets_;
    other.uninit();
}

template<class T, int R>
shared_shape<T,R>::~shared_shape()
{
    decref();
}

template<class T, int R>
void shared_shape<T,R>::copy_before_write()
{
    if (R > 1 and refs_ and *refs_ > 1) 
        *this = this->copy();
}

template<class T, int R>
struct _data_from_ptrs_noffsets_ndataoffsets {
    static
    T* call(typename PointerArray<T,R>::type ptrs,
            typename shared_shape<T,R>::size_type noffsets,
            typename shared_shape<T,R>::size_type ndataoffsets)
    {
        return (T*)(ptrs[noffsets - ndataoffsets]);
    }
};

template<class T>
struct _data_from_ptrs_noffsets_ndataoffsets<T,1> {
    static T* call (typename PointerArray<T,1>::type ptrs,
                    typename shared_shape<T,1>::size_type noffsets,
                    typename shared_shape<T,1>::size_type ndataoffsets)
    {
        return (T*)(ptrs);
    }
};

template<class T, int R>
T* shared_shape<T,R>::data() const
{
    return _data_from_ptrs_noffsets_ndataoffsets<T,R>::call(ptrs_,noffsets_,ndataoffsets_);
}

template<class T, int R>
typename shared_shape<T,R>::size_type shared_shape<T,R>::size() const
{
    return ndataoffsets_ * extent_[R-1];
}


template<class T, int R>
typename shared_shape<T,R>::size_type shared_shape<T,R>::extent(int i) const
{
    return extent_[i];
}

template<class T, int R> 
const std::array<typename shared_shape<T,R>::size_type,R>& shared_shape<T,R>::extent() const {
    return extent_;
}

template<class T, int R>
void shared_shape<T,R>::relocate(T* newdata)
{
    if (R==1) {
        ptrs_ = reinterpret_cast<ptrs_type>(newdata);
    } else if (R>1) {
        std::ptrdiff_t shift = newdata - data();
        if (shift != 0) {
            copy_before_write();
            for (size_type i = noffsets_ - ndataoffsets_; i < noffsets_; i++)
                orig_[i] = reinterpret_cast<void**>((T*)(orig_[i]) + shift);
        }
    }
}

template<class T, int R>
shared_shape<T,R> shared_shape<T,R>::copy() const
{
    shared_shape copy_of_this;
    copy_of_this.extent_ = extent_;
    copy_of_this.noffsets_ = noffsets_;
    copy_of_this.ndataoffsets_ = ndataoffsets_;
    if (R>1) {
        copy_of_this.refs_ = new std::atomic<int>(1);
        copy_of_this.orig_ = new void**[noffsets_];
        std::copy((void***)ptrs_, (void***)(ptrs_) + noffsets_, copy_of_this.orig_);
        std::ptrdiff_t shift = copy_of_this.orig_ - (void***)(ptrs_);
        for (size_type i = 0; i < noffsets_ - ndataoffsets_; i++)
            copy_of_this.orig_[i] += shift;
        copy_of_this.ptrs_ = reinterpret_cast<ptrs_type>(copy_of_this.orig_);
    } else {        
        copy_of_this.ptrs_ = ptrs_;
    }
    return copy_of_this;
}

template<class T, int R>
shared_shape<T,R-1> shared_shape<T,R>::at(size_type index) const
{
    if (R < 1 or index < 0 or index >= extent_[0])
        throw std::out_of_range("shared_shape::at");
    shared_shape<T,R-1> result;
    if (R>1) {
        for (size_type i = 0; i < R-1; ++i)
            result.extent_[i] = extent_[i+1];
        result.ptrs_         = ptrs_[index];
        result.refs_         = refs_;
        result.orig_         = orig_;
        result.noffsets_     = noffsets_/extent_[0] - 1;
        result.ndataoffsets_ = ndataoffsets_/extent_[0];
        incref();
    }
    return result;
}

template<class T, int R>
void shared_shape<T,R>::reshape(const std::array<size_type,R>&extent)
{
    if (extent != extent_) {
        *this = shared_shape<T,R>(extent, data());
    }
}

template<class T, int R>
typename shared_shape<T,R>::ptrs_type shared_shape<T,R>::ptrs() const
{
    return ptrs_;
}

template<class T, int R>
void shared_shape<T,R>::uninit()
{
    ptrs_         = nullptr;
    orig_         = nullptr;
    refs_         = nullptr;
    noffsets_     = 0;
    ndataoffsets_ = 0;
    extent_.fill(0);
}

template<class T, int R>
void shared_shape<T,R>::incref() const
{
    if (refs_)
        (*refs_)++;
}

template<class T, int R>
void shared_shape<T,R>::decref()
{
    if (refs_) {
        if (--(*refs_) == 0) {
            if (R>1)
                delete[] orig_;
            delete refs_;
            uninit();
        }
    }
}

/***************************************************************************/

}

#endif


namespace ra {

typedef ssize_t size_type;
enum class RESIZE { NO, ALLOWED };

}


namespace ra { template<typename T,int R> class rarray;    }
namespace ra { template<typename T>       class CommaOp;   }


#define ExOp class

namespace ra { 
template<typename T, int R, ExOp AOP, typename A1, typename A2, typename A3> class Expr;
}
   
namespace ra {

    
template<typename T,int R> 
class rarray {
    
  public: 
    typedef ssize_t                                   difference_type;  // difference type for indices
    typedef ssize_t                                   size_type;        // type of indices
    typedef T*                                        iterator;         // iterator type
    typedef const T*                                  const_iterator;   // iterator type for constant access
    typedef typename PointerArray<T,R>::type          parray_t;         // shorthand for T*const*const*...
    typedef typename PointerArray<T,R>::noconst_type  noconst_parray_t; // shorthand for T***...
    RA_INLINE_ rarray();                                                                                                                  // constructor leaving rarray undefined 
    RA_INLINE_ explicit rarray(size_type n0);                                                                                             // constructor creating its own buffer for R=1
    RA_INLINE_ rarray(size_type n0, size_type n1);                                                                                        // constructor creating its own buffer for R=2
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2);                                                                                                                  // R=3
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2, size_type n3);                                                                                                    // R=4
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4);                                                                                      // R=5
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5);                                                                        // R=6
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6);                                                          // R=7
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6, size_type n7);                                            // R=8
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6, size_type n7, size_type n8);                              // R=9 
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6, size_type n7, size_type n8, size_type n9);                // R=10
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6, size_type n7, size_type n8, size_type n9, size_type n10); // R=11
    RA_INLINE_ rarray(const size_type* extent);                                                                                                                                   // R>11
    RA_INLINE_ rarray(T* buffer, size_type n0);                                                                                           // constructor from an existing buffer for R=1
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1);                                                                                                                     // R=2
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1, size_type n2);                                                                                                       // R=3
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1, size_type n2, size_type n3);                                                                                         // R=4
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1, size_type n2, size_type n3, size_type n4);                                                                           // R=5
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5);                                                             // R=6
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6);                                               // R=7 
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6, size_type n7);                                 // R=8
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6, size_type n7, size_type n8);                   // R=9
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6, size_type n7, size_type n8, size_type n9);     // R=10
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6, size_type n7, size_type n8, size_type n9, size_type n10);// R=11
    RA_INLINE_ rarray(T* buffer, const size_type* extent);                                                                                                                        // R>11
    RA_INLINEF rarray(const rarray<T,R> &a);                                      // copy constructor
    RA_INLINE_ rarray<T,R>& operator=(const rarray<T,R> &a);                      // array assignment operator
    rarray(rarray<T,R>&& x);                                                      // move constructor
    rarray<T,R>& operator=(rarray<T,R>&& x);                                      // move assignment operator
    RA_INLINE_ CommaOp<T> operator=(const T& e);                                  // Comma separated element assignment
    RA_INLINEF ~rarray();                                                         // destructor
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF explicit   rarray (const Expr<T,R,AOP,A1,A2,A3>& e);
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF rarray& operator= (const Expr<T,R,AOP,A1,A2,A3>& e);
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF rarray& operator+=(const Expr<T,R,AOP,A1,A2,A3>& e);
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF rarray& operator-=(const Expr<T,R,AOP,A1,A2,A3>& e);
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF rarray& operator*=(const Expr<T,R,AOP,A1,A2,A3>& e);
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF rarray& operator/=(const Expr<T,R,AOP,A1,A2,A3>& e);
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF rarray& operator%=(const Expr<T,R,AOP,A1,A2,A3>& e);
    RA_INLINEF void clear();                                                      // clean up routine, make undefined
    RA_INLINE_ void reshape(size_type n0, RESIZE resize_allowed=RESIZE::NO);                                                                            // reshape shallow copy keeping the underlying data for R=1
    RA_INLINE_ void reshape(size_type n0, size_type n1, RESIZE resize_allowed=RESIZE::NO);                                                                                                                   // R=2
    RA_INLINE_ void reshape(size_type n0, size_type n1, size_type n2, RESIZE resize_allowed=RESIZE::NO);                                                                                                     // R=3
    RA_INLINE_ void reshape(size_type n0, size_type n1, size_type n2, size_type n3, RESIZE resize_allowed=RESIZE::NO);                                                                                       // R=4
    RA_INLINE_ void reshape(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, RESIZE resize_allowed=RESIZE::NO);                                                                         // R=5
    RA_INLINE_ void reshape(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, RESIZE resize_allowed=RESIZE::NO);                                                           // R=6
    RA_INLINE_ void reshape(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6, RESIZE resize_allowed=RESIZE::NO);                                             // R=7
    RA_INLINE_ void reshape(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6, size_type n7, RESIZE resize_allowed=RESIZE::NO);                               // R=8
    RA_INLINE_ void reshape(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6, size_type n7, size_type n8, RESIZE resize_allowed=RESIZE::NO);                 // R=9
    RA_INLINE_ void reshape(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6, size_type n7, size_type n8, size_type n9, RESIZE resize_allowed=RESIZE::NO);   // R=10
    RA_INLINE_ void reshape(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6, size_type n7, size_type n8, size_type n9, size_type n10, RESIZE resize_allowed=RESIZE::NO); // R=11
    RA_INLINE_ void reshape(const size_type* extent, RESIZE resize_allowed=RESIZE::NO);                                                                                                                      // R>11
    RA_INLINE_ bool                is_clear()           const;                       // check if undefined
    RA_INLINE_ rarray<T,R>         copy()               const;                       // return a copy
    RA_INLINE_ size_type           extent(int i)        const;                       // retrieve array size in dimension i
    RA_INLINE_ const size_type*    shape()              const;                       // retrieve array sizes in all dimensions
    RA_INLINE_ size_type           size()               const;                       // retrieve the total number of elements  
    RA_INLINE_ T*                  data();                                           // return a T* to the internal data
    RA_INLINE_ const T*            data()               const;                       // return a T* to the internal data
    RA_INLINE_ parray_t            ptr_array()          const;                       // return a T*const*.. acting similarly to this rarray when using []:
    RA_INLINE_ noconst_parray_t    noconst_ptr_array()  const;                       // return a T**.. acting similarly to this rarray when using []:    
    RA_INLINE_ rarray<const T,R>&  const_ref()          const;                       // create a reference to this that treats elements as constant:
    RA_INLINE_ void                fill(const T& value);                             // fill with uniform value
    RA_INLINE_ iterator            begin();                                          // start of the content
    RA_INLINE_ const_iterator      begin()              const;                       // start of the content, when *this is constant
    RA_INLINE_ const_iterator      cbegin()             const;                       // start of the content, when *this can be constant and you need to be explicit
    RA_INLINE_ iterator            end();                                            // end of the content
    RA_INLINE_ const_iterator      end()                const;                       // end of the content, when *this is constant
    RA_INLINE_ const_iterator      cend()               const;                       // end of the content, when *this is constant and you need to be explicit about that
    RA_INLINE_ size_type           index(const T& a, int i) const;                   // if a an element in the array, get index i of that element
    RA_INLINE_ size_type           index(const iterator& iter, int i) const;         // if i points at an element in the array, get index i of that element
    RA_INLINE_ std::array<size_type,R> index(const T& a) const;           // if a an element in the array, get the indices of that element
    RA_INLINE_ std::array<size_type,R> index(const iterator& i) const;// if i points at an element in the array, get the indices of that element
    RA_INLINE_ size_type*          index(const T& a, size_type* ind) const;           // if a an element in the array, get the indices of that element
    RA_INLINE_ size_type*          index(const iterator& i, size_type* ind) const;          // if i points at an element in the array, get the indices of that element
    RA_INLINEF rarray<T,R-1> at(size_type i);
    RA_INLINEF rarray<T,R-1> at(size_type i) const;
    RA_INLINEF operator typename PointerArray<T,R>::type ();  // makes a[..][..] work, as well as automatic conversion
    RA_INLINEF operator typename PointerArray<const T,R>::type () const; 
    RA_INLINEF const T& leval(size_type i) const;
    rarray(shared_buffer<T>&& buffer, shared_shape<T,R>&& shape) :
        buffer_(std::forward<shared_buffer<T>>(buffer)),
        shape_(std::forward<shared_shape<T,R>>(shape))
    {}
  private:
    shared_buffer<T>  buffer_;
    shared_shape<T,R> shape_;
        
}; // end definition rarray<T,R>

template<typename T> 
class rarray<T,0> {
    
  public: 
    typedef ssize_t                                   difference_type;  // difference type for indices
    typedef ssize_t                                   size_type;        // type of indices
    typedef T*                                        iterator;         // iterator type
    typedef const T*                                  const_iterator;   // iterator type for constant access
    typedef typename PointerArray<T,0>::type          parray_t;         // shorthand for T*const*const*...
    typedef typename PointerArray<T,0>::noconst_type  noconst_parray_t; // shorthand for T***...
    RA_INLINE_ rarray() {};                                                                                                                  // constructor leaving rarray undefined 
    RA_INLINE_ rarray(const size_type* extent);                                                                                                                                   // R>11
    RA_INLINE_ rarray(T* buffer, const size_type* extent);                                                                                                                        // R>11
    RA_INLINEF rarray(const rarray<T,0> &a);                                      // copy constructor
    RA_INLINE_ rarray<T,0>& operator=(const rarray<T,0> &a);                      // array assignment operator
    rarray(rarray<T,0>&& x);                                                      // move constructor
    rarray<T,0>& operator=(rarray<T,0>&& x);                                      // move assignment operator
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF explicit   rarray (const Expr<T,0,AOP,A1,A2,A3>& e);
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF rarray& operator= (const Expr<T,0,AOP,A1,A2,A3>& e);
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF rarray& operator+=(const Expr<T,0,AOP,A1,A2,A3>& e);
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF rarray& operator-=(const Expr<T,0,AOP,A1,A2,A3>& e);
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF rarray& operator*=(const Expr<T,0,AOP,A1,A2,A3>& e);
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF rarray& operator/=(const Expr<T,0,AOP,A1,A2,A3>& e);
    template<ExOp AOP, typename A1, typename A2, typename A3> RA_INLINEF rarray& operator%=(const Expr<T,0,AOP,A1,A2,A3>& e);
    RA_INLINEF void clear();                                                      // clean up routine, make undefined
    RA_INLINE_ bool                is_clear()           const {  return buffer_.cbegin() == nullptr; }                      // check if undefined
    RA_INLINE_ rarray<T,0>         copy()               const;                       // return a copy
    RA_INLINE_ size_type           extent(int i)        const { return is_clear()?0:1; }   // retrieve array size in dimension i
    RA_INLINE_ const size_type*       shape()              const;                       // retrieve array sizes in all dimensions
    RA_INLINE_ size_type           size()               const;                       // retrieve the total number of elements  
    RA_INLINE_ T*                  data();                                           // return a T* to the internal data
    RA_INLINE_ const T*            data()               const { return buffer_.begin(); }                   // return a T* to the internal data
    RA_INLINE_ parray_t            ptr_array()          const;                       // return a T*const*.. acting similarly to this rarray when using []:
    RA_INLINE_ noconst_parray_t    noconst_ptr_array()  const;                       // return a T**.. acting similarly to this rarray when using []:    
    RA_INLINE_ rarray<const T,0>&  const_ref()          const;                       // create a reference to this that treats elements as constant:
    RA_INLINE_ void                fill(const T& value);                             // fill with uniform value
    RA_INLINE_ iterator            begin();                                          // start of the content
    RA_INLINE_ const_iterator      begin()              const;                       // start of the content, when *this is constant
    RA_INLINE_ const_iterator      cbegin()             const;                       // start of the content, when *this can be constant and you need to be explicit
    RA_INLINE_ iterator            end();                                            // end of the content
    RA_INLINE_ const_iterator      end()                const;                       // end of the content, when *this is constant
    RA_INLINE_ const_iterator      cend()               const;                       // end of the content, when *this is constant and you need to be explicit about that
    RA_INLINE_ size_type           index(const T& a, size_type i) const;             // if a an element in the array, get index i of that element
    RA_INLINE_ size_type           index(const iterator& iter, size_type i);         // if i points at an element in the array, get index i of that element
    RA_INLINE_ size_type           index(const const_iterator& iter, size_type i) const;// if i points at an element in the array, get index i of that element
    RA_INLINE_ size_type*          index(const T& a, size_type* ind) const;           // if a an element in the array, get the indices of that element
    RA_INLINE_ size_type*          index(const iterator& i, size_type* ind) const;          // if i points at an element in the array, get the indices of that element
    RA_INLINE_ size_type*          index(const const_iterator& i, size_type* ind) const;// if i points at an element in the array, get the indices of that element
    RA_INLINEF const T& leval(size_type i) const;
    rarray(shared_buffer<T>&& buffer, const shared_shape<T,0>& shape) :
        buffer_(std::forward<shared_buffer<T>>(buffer))
    {}
    ~rarray() {}
    RA_INLINEF operator T& () { return buffer_[0]; }
    RA_INLINEF operator const T& () const { return buffer_[0]; }
    RA_INLINEF T& at(size_type i) { return buffer_[0]; }
    RA_INLINEF const T& at(size_type i) const { return buffer_[0]; }
    RA_INLINE_ CommaOp<T> operator=(const T& e)                                 // Comma separated element assignment
    {
        T* first = &(buffer_[0]);
        *first = e;
        ra::CommaOp<T> co(first+1, first); 
        return co;
    }
  private:
    shared_buffer<T>  buffer_;
}; // end definition rarray<T,0>

template<typename T>
class CommaOp {
  public:
    RA_INLINEF CommaOp& operator,(const T& e);                         // puts the next number into the array.
  private:
    RA_INLINEF CommaOp(T* ptr, T* last); 
    T *ptr_;                                                           // points to next element to be filled
    T * const last_;                                                   // points to last element
    template<typename,int> friend class rarray;
    template<typename,int> friend class subrarray;
};

template<typename A>                                                              RA_INLINE_ int extent_given_byte_size(A a[], int i, int byte_size);                               //for 1d array
template<typename A,int Z>                                                        RA_INLINE_ int extent_given_byte_size(A a[][Z], int i, int byte_size);                            //for 2d array
template<typename A,int Y,int Z>                                                  RA_INLINE_ int extent_given_byte_size(A a[][Y][Z], int i, int byte_size);                         //for 3d array
template<typename A,int X,int Y,int Z>                                            RA_INLINE_ int extent_given_byte_size(A a[][X][Y][Z], int i, int byte_size);                      //for 4d array
template<typename A,int W,int X,int Y,int Z>                                      RA_INLINE_ int extent_given_byte_size(A a[][W][X][Y][Z], int i, int byte_size);                   //for 5d array
template<typename A,int V,int W,int X,int Y,int Z>                                RA_INLINE_ int extent_given_byte_size(A a[][V][W][X][Y][Z], int i, int byte_size);                //for 6d array
template<typename A,int U,int V,int W,int X,int Y,int Z>                          RA_INLINE_ int extent_given_byte_size(A a[][U][V][W][X][Y][Z], int i, int byte_size);             //for 7d array
template<typename A,int T,int U,int V,int W,int X,int Y,int Z>                    RA_INLINE_ int extent_given_byte_size(A a[][T][U][V][W][X][Y][Z], int i, int byte_size);          //for 8d array
template<typename A,int S,int T,int U,int V,int W,int X,int Y,int Z>              RA_INLINE_ int extent_given_byte_size(A a[][S][T][U][V][W][X][Y][Z], int i, int byte_size);       //for 9d array
template<typename A,int R,int S,int T,int U,int V,int W,int X,int Y,int Z>        RA_INLINE_ int extent_given_byte_size(A a[][R][S][T][U][V][W][X][Y][Z], int i, int byte_size);    //for 10d array
template<typename A,int Q,int R,int S,int T,int U,int V,int W,int X,int Y,int Z>  RA_INLINE_ int extent_given_byte_size(A a[][Q][R][S][T][U][V][W][X][Y][Z], int i, int byte_size); //for 11d array
template<typename A,int R>                                                        RA_INLINE_ int extent_given_byte_size(const rarray<A,R>& a, int i, int byte_size); // use rarray's extent function
template<typename A>                                                              RA_INLINE_ rarray<A,1>  make_rarray_given_byte_size(A a[], int byte_size);                              //for 1d array
template<typename A,int Z>                                                        RA_INLINE_ rarray<A,2>  make_rarray_given_byte_size(A a[][Z], int byte_size);                           //for 2d array
template<typename A,int Y,int Z>                                                  RA_INLINE_ rarray<A,3>  make_rarray_given_byte_size(A a[][Y][Z], int byte_size);                        //for 3d array
template<typename A,int X,int Y,int Z>                                            RA_INLINE_ rarray<A,4>  make_rarray_given_byte_size(A a[][X][Y][Z], int byte_size);                     //for 4d array
template<typename A,int W,int X,int Y,int Z>                                      RA_INLINE_ rarray<A,5>  make_rarray_given_byte_size(A a[][W][X][Y][Z], int byte_size);                  //for 5d array
template<typename A,int V,int W,int X,int Y,int Z>                                RA_INLINE_ rarray<A,6>  make_rarray_given_byte_size(A a[][V][W][X][Y][Z], int byte_size);               //for 6d array
template<typename A,int U,int V,int W,int X,int Y,int Z>                          RA_INLINE_ rarray<A,7>  make_rarray_given_byte_size(A a[][U][V][W][X][Y][Z], int byte_size);            //for 7d array
template<typename A,int T,int U,int V,int W,int X,int Y,int Z>                    RA_INLINE_ rarray<A,8>  make_rarray_given_byte_size(A a[][T][U][V][W][X][Y][Z], int byte_size);         //for 8d array
template<typename A,int S,int T,int U,int V,int W,int X,int Y,int Z>              RA_INLINE_ rarray<A,9>  make_rarray_given_byte_size(A a[][S][T][U][V][W][X][Y][Z], int byte_size);      //for 9d array
template<typename A,int R,int S,int T,int U,int V,int W,int X,int Y,int Z>        RA_INLINE_ rarray<A,10> make_rarray_given_byte_size(A a[][R][S][T][U][V][W][X][Y][Z], int byte_size);   //for 10d array
template<typename A,int Q,int R,int S,int T,int U,int V,int W,int X,int Y,int Z>  RA_INLINE_ rarray<A,11> make_rarray_given_byte_size(A a[][Q][R][S][T][U][V][W][X][Y][Z], int byte_size);//for 11d array
template<typename A,int R>                                                        RA_INLINE_ rarray<A,R> make_rarray_given_byte_size(rarray<A,R> a, int byte_size); // trivial action for rarray

template<typename S>
rarray<S,1> linspace(S x1, S x2, int n=0, bool end_incl=true)
{
    if (n==0) {
        if (x2>x1)
            n = x2 - x1 + end_incl;
        else
            n = x1 - x2 + end_incl;
    }
    rarray<S,1> x(n);
    for (int i = 0; i < n; i++)
        x[i] = x1 + ((x2-x1)*(long long int)(i))/(n-end_incl);
    if (end_incl)
        x[n-1] = x2;
    return x;
}

template<class T>
class Xrange {
  private:
    struct const_iterator {
        typedef std::ptrdiff_t difference_type;
        typedef T value_type;
        typedef T* pointer;
        typedef T& reference;
        typedef std::input_iterator_tag iterator_category;
        const_iterator(T i, T di, T b): i_(i), di_(di), b_(b) {}
        bool operator!=(const const_iterator& other) const {
            return i_ != other.i_;
        }
        const const_iterator& operator++() {
            i_+=di_;
            if (di_>0 && i_ >= b_)
               i_ = b_;
            if (di_<0 && i_ <= b_)
               i_ = b_;
            return *this;
        }
        const T& operator*() {
            return i_;
        }
        T i_, di_, b_;
    };
    T a_, b_, d_;
  public:
    Xrange(T a, T b, T d)
      : a_(a), b_(a + size_t(std::ceil((b-a)/double(d)))*d), d_(d) 
    {}
    const_iterator begin() const {
        return const_iterator(a_, d_, b_);
    }
    const_iterator end() const {
        return const_iterator(b_, d_, b_);
    }
};

template<class T>
Xrange<T> xrange(T end) {
    return Xrange<T>((T)0, end, (T)1);
}

template<class S, class T>
Xrange<T> xrange(S begin, T end) {
    return Xrange<T>((T)begin, end, (T)1);
}

template<class S, class T, class U>
Xrange<T> xrange(S begin, T end, U step) {
    return Xrange<T>((T)begin, end, (T)step);
}


} // end namespace ra

#include <utility>

#include <algorithm>

    

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray()
  : buffer_(),
    shape_()
{}
    
template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(ra::size_type n0)
  : buffer_(n0),
    shape_({n0}, buffer_.begin())
{
    static_assert( R==1, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(size_type n0,size_type n1)
  : buffer_(n0*n1),
    shape_({n0,n1}, buffer_.begin())
{
    static_assert( R==2, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(size_type n0,size_type n1,size_type n2)
  : buffer_(n0*n1*n2),
    shape_({n0,n1,n2}, buffer_.begin())
{
    static_assert( R==3, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(size_type n0,size_type n1,size_type n2,size_type n3)
  : buffer_(n0*n1*n2*n3),
    shape_({n0,n1,n2,n3}, buffer_.begin())
{
    static_assert( R==4, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4)
  : buffer_(n0*n1*n2*n3*n4),
    shape_({n0,n1,n2,n3,n4}, buffer_.begin())
{
    static_assert( R==5, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5)
  : buffer_(n0*n1*n2*n3*n4*n5),
    shape_({n0,n1,n2,n3,n4,n5}, buffer_.begin())
{
    static_assert( R==6, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6)
  : buffer_(n0*n1*n2*n3*n4*n5*n6),
    shape_({n0,n1,n2,n3,n4,n5,n6}, buffer_.begin())
{
    static_assert( R==7, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6,size_type n7)
  : buffer_(n0*n1*n2*n3*n4*n5*n6*n7),
    shape_({n0,n1,n2,n3,n4,n5,n6,n7}, buffer_.begin())
{
    static_assert( R==8, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6,size_type n7,size_type n8)
  : buffer_(n0*n1*n2*n3*n4*n5*n6*n7*n8),
    shape_({n0,n1,n2,n3,n4,n5,n6,n7,n8}, buffer_.begin())
{
    static_assert( R==9, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6,size_type n7,size_type n8,size_type n9)
  : buffer_(n0*n1*n2*n3*n4*n5*n6*n7*n8*n9),
    shape_({n0,n1,n2,n3,n4,n5,n6,n7,n8,n9}, buffer_.begin())
{
    static_assert( R==10, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6,size_type n7,size_type n8,size_type n9,size_type n10)
  : buffer_(n0*n1*n2*n3*n4*n5*n6*n7*n8*n9*n10),
    shape_({n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10}, buffer_.begin())
{
    static_assert( R==11, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(T* buffer, size_type n0)
  : buffer_(n0,  buffer),
    shape_({n0}, buffer)
{
    static_assert( R==1, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(T* buffer, size_type n0,size_type n1)
  : buffer_(n0*n1,  buffer),
    shape_({n0,n1}, buffer)
{
    static_assert( R==2, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(T* buffer, size_type n0,size_type n1,size_type n2)
  : buffer_(n0*n1*n2,  buffer),
    shape_({n0,n1,n2}, buffer)
{
    static_assert( R==3, "Incorrect number of dimensions in rarray constructor.");
}


template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(T* buffer, size_type n0,size_type n1,size_type n2,size_type n3)
  : buffer_(n0*n1*n2*n3,  buffer),
    shape_({n0,n1,n2,n3}, buffer)
{
    static_assert( R==4, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(T* buffer, size_type n0,size_type n1,size_type n2,size_type n3,size_type n4)
  : buffer_(n0*n1*n2*n3*n4,  buffer),
    shape_({n0,n1,n2,n3,n4}, buffer)
{
    static_assert( R==5, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(T* buffer, size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5)
  : buffer_(n0*n1*n2*n3*n4*n5,  buffer),
    shape_({n0,n1,n2,n3,n4,n5}, buffer)
{
    static_assert( R==6, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(T* buffer, size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6)
  : buffer_(n0*n1*n2*n3*n4*n5*n6,  buffer),
    shape_({n0,n1,n2,n3,n4,n5,n6}, buffer)
{
    static_assert( R==7, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(T* buffer, size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6,size_type n7)
  : buffer_(n0*n1*n2*n3*n4*n5*n6*n7,  buffer),
    shape_({n0,n1,n2,n3,n4,n5,n6,n7}, buffer)
{
    static_assert( R==8, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(T* buffer, size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6,size_type n7,size_type n8)
  : buffer_(n0*n1*n2*n3*n4*n5*n6*n7*n8,  buffer),
    shape_({n0,n1,n2,n3,n4,n5,n6,n7,n8}, buffer)
{
    static_assert( R==9, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(T* buffer, size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6,size_type n7,size_type n8,size_type n9)
  : buffer_(n0*n1*n2*n3*n4*n5*n6*n7*n8*n9,  buffer),
    shape_({n0,n1,n2,n3,n4,n5,n6,n7,n8,n9}, buffer)
{
    static_assert( R==10, "Incorrect number of dimensions in rarray constructor.");
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(T* buffer, size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6,size_type n7,size_type n8,size_type n9,size_type n10)
  : buffer_(n0*n1*n2*n3*n4*n5*n6*n7*n8*n9*n10,  buffer),
    shape_({n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10}, buffer)
{
    static_assert( R==11, "Incorrect number of dimensions in rarray constructor.");
}

namespace ra {
RA_INLINE_ size_type mul(const size_type * x, int n) {
    size_type result = 1;
    for (int i=0;i<n;i++)
        result *= x[i];
    return result;
}
}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(const size_type* extent)
  : buffer_(mul(extent,R)),
    shape_((const std::array<size_type,R>&)(*extent), buffer_.begin())
{}

template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::rarray(T* buffer, const size_type* extent)
  : buffer_(mul(extent,R), buffer),
    shape_((const std::array<size_type,R>&)(*extent), buffer)
{}

template<typename T, int R> RA_INLINEF
ra::rarray<T,R>::rarray(const rarray<T,R> &a)
  : buffer_(a.buffer_),
    shape_(a.shape_)
{
}

template<typename T, int R> RA_INLINEF
ra::rarray<T,R>::rarray(rarray<T,R>&& x)
  : buffer_(std::move(x.buffer_)),
    shape_(std::move(x.shape_))
{
}


template<typename T, int R> RA_INLINE_ 
void ra::rarray<T,R>::clear()
{
    shape_ = shared_shape<T,R>();
    buffer_ = shared_buffer<T>();
}


template<typename T, int R> RA_INLINE_ 
void ra::rarray<T,R>::reshape(size_type n0, RESIZE resize_allowed)
{
    static_assert( R==1, "Incorrect number of dimensions in rarray::reshape method.");
    
    if (size() == n0
        or (resize_allowed == RESIZE::ALLOWED and size() >= n0)) {
        
        shape_ = shared_shape<T,R>({n0}, buffer_.begin());

    } else {

        throw std::out_of_range(std::string("Incompatible dimensions in function ") + std::string(__PRETTY_FUNCTION__));
        
    }
}

template<typename T, int R> RA_INLINE_ 
void ra::rarray<T,R>::reshape(size_type n0, size_type n1, RESIZE resize_allowed)
{
    static_assert( R==2, "Incorrect number of dimensions in rarray::reshape method.");
    
    if (size() == n0*n1
        or (resize_allowed == RESIZE::ALLOWED and size() >= n0*n1)) {
        
        shape_ = shared_shape<T,R>({n0,n1}, buffer_.begin());
        
    } else {
        
        throw std::out_of_range(std::string("Incompatible dimensions in function ") + std::string(__PRETTY_FUNCTION__));
        
    }
}

template<typename T, int R> RA_INLINE_ 
void ra::rarray<T,R>::reshape(size_type n0, size_type n1, size_type n2, RESIZE resize_allowed)
{
    static_assert( R==3, "Incorrect number of dimensions in rarray::reshape method.");
    
    if (size() == n0*n1*n2
        or (resize_allowed == RESIZE::ALLOWED and size() >= n0*n1*n2)) {
        
        shape_ = shared_shape<T,R>({n0,n1,n2}, buffer_.begin());
        
    } else {
            
        throw std::out_of_range(std::string("Incompatible dimensions in function ") + std::string(__PRETTY_FUNCTION__));
        
    }
}


template<typename T, int R> RA_INLINE_ 
void ra::rarray<T,R>::reshape(size_type n0,size_type n1,size_type n2,size_type n3, RESIZE resize_allowed)
{
    static_assert( R==4, "Incorrect number of dimensions in rarray::reshape method.");

    if (size() == n0*n1*n2*n3
        or (resize_allowed == RESIZE::ALLOWED and size() >= n0*n1*n2*n3)) {
        
        shape_ = shared_shape<T,R>({n0,n1,n2,n3}, buffer_.begin());
        
    } else {
        
        throw std::out_of_range(std::string("Incompatible dimensions in function ") + std::string(__PRETTY_FUNCTION__));
        
    }
}

template<typename T, int R> RA_INLINE_ 
void ra::rarray<T,R>::reshape(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4, RESIZE resize_allowed)
{
    static_assert( R==5, "Incorrect number of dimensions in rarray::reshape method.");
    
    if (size() == n0*n1*n2*n3*n4
        or (resize_allowed == RESIZE::ALLOWED and size() >= n0*n1*n2*n3*n4)) {
        
        shape_ = shared_shape<T,R>({n0,n1,n2,n3,n4}, buffer_.begin());
        
    } else {
        
        throw std::out_of_range(std::string("Incompatible dimensions in function ") + std::string(__PRETTY_FUNCTION__));
        
    }
}

template<typename T, int R> RA_INLINE_ 
void ra::rarray<T,R>::reshape(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5, RESIZE resize_allowed)
{
    static_assert( R==6, "Incorrect number of dimensions in rarray::reshape method.");

    if (size() == n0*n1*n2*n3*n4*n5
        or (resize_allowed == RESIZE::ALLOWED and size() >= n0*n1*n2*n3*n4*n5)) {
        
        shape_ = shared_shape<T,R>({n0,n1,n2,n3,n4,n5}, buffer_.begin());
        
    } else {
        
        throw std::out_of_range(std::string("Incompatible dimensions in function ") + std::string(__PRETTY_FUNCTION__));
        
    }
}

template<typename T, int R> RA_INLINE_ 
void ra::rarray<T,R>::reshape(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6, RESIZE resize_allowed)
{
    static_assert( R==7, "Incorrect number of dimensions in rarray::reshape method.");
        
    if (size() == n0*n1*n2*n3*n4*n5*n6
        or (resize_allowed == RESIZE::ALLOWED and size() >= n0*n1*n2*n3*n4*n5*n6)) {
        
        shape_ = shared_shape<T,R>({n0,n1,n2,n3,n4,n5,n6}, buffer_.begin());
        
    } else {
        
        throw std::out_of_range(std::string("Incompatible dimensions in function ") + std::string(__PRETTY_FUNCTION__));
        
    }
}

template<typename T, int R> RA_INLINE_ 
void ra::rarray<T,R>::reshape(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6,size_type n7, RESIZE resize_allowed)
{
    static_assert( R==8, "Incorrect number of dimensions in rarray::reshape method.");
        
    if (size() == n0*n1*n2*n3*n4*n5*n6*n7
        or (resize_allowed == RESIZE::ALLOWED and size() >= n0*n1*n2*n3*n4*n5*n6*n7)) {
        
        shape_ = shared_shape<T,R>({n0,n1,n2,n3,n4,n5,n6,n7}, buffer_.begin());
        
    } else {
        
        throw std::out_of_range(std::string("Incompatible dimensions in function ") + std::string(__PRETTY_FUNCTION__));
        
    }
}

template<typename T, int R> RA_INLINE_ 
void ra::rarray<T,R>::reshape(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6,size_type n7,size_type n8, RESIZE resize_allowed)
{
    static_assert( R==9, "Incorrect number of dimensions in rarray::reshape method.");
    
    if (size() == n0*n1*n2*n3*n4*n5*n6*n7*n8
        or (resize_allowed == RESIZE::ALLOWED and size() >= n0*n1*n2*n3*n4*n5*n6*n7*n8)) {
        
        shape_ = shared_shape<T,R>({n0,n1,n2,n3,n4,n5,n6,n7,n8}, buffer_.begin());
        
    } else {
        
        throw std::out_of_range(std::string("Incompatible dimensions in function ") + std::string(__PRETTY_FUNCTION__));
        
    }
}

template<typename T, int R> RA_INLINE_ 
void ra::rarray<T,R>::reshape(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6,size_type n7,size_type n8,size_type n9, RESIZE resize_allowed)
{
    static_assert( R==10, "Incorrect number of dimensions in rarray::reshape method.");
    
    if (size() == n0*n1*n2*n3*n4*n5*n6*n7*n8*n9
        or (resize_allowed == RESIZE::ALLOWED and size() >= n0*n1*n2*n3*n4*n5*n6*n7*n8*n9)) {
        
        shape_ = shared_shape<T,R>({n0,n1,n2,n3,n4,n5,n6,n7,n8,n9}, buffer_.begin());
        
    } else {
        
        throw std::out_of_range(std::string("Incompatible dimensions in function ") + std::string(__PRETTY_FUNCTION__));
        
    }
}

template<typename T, int R> RA_INLINE_ 
void ra::rarray<T,R>::reshape(size_type n0,size_type n1,size_type n2,size_type n3,size_type n4,size_type n5,size_type n6,size_type n7,size_type n8,size_type n9,size_type n10, RESIZE resize_allowed)
{
    static_assert( R==11, "Incorrect number of dimensions in rarray::reshape method.");
    
    if (size() == n0*n1*n2*n3*n4*n5*n6*n7*n8*n9*n10
        or (resize_allowed == RESIZE::ALLOWED and size() >= n0*n1*n2*n3*n4*n5*n6*n7*n8*n9*n10)) {
        
        shape_ = shared_shape<T,R>({n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10}, buffer_.begin());

    } else {
        
        throw std::out_of_range(std::string("Incompatible dimensions in function ") + std::string(__PRETTY_FUNCTION__));

    }
}

template<typename T, int R> RA_INLINE_ 
void ra::rarray<T,R>::reshape(const size_type* newshape, RESIZE resize_allowed)
{
    if (size() == mul(newshape,R)
        or (resize_allowed == RESIZE::ALLOWED and size() >= mul(newshape,R))) {
        
        shape_ = shared_shape<T,R>((const std::array<size_type,R>&)(*newshape), buffer_.begin());
        
    } else {
        
        throw std::out_of_range(std::string("Incompatible dimensions in function ") + std::string(__PRETTY_FUNCTION__));
        
    }
}



template<typename T, int R> RA_INLINE_ 
ra::rarray<T,R>::~rarray()
{}

template<typename T, int R> RA_INLINE_
ra::CommaOp<T> ra::rarray<T,R>::operator=(const T& e)
{
    RA_CHECKORSAY(parray_!=nullptr, "assignment to unsized array");
    RA_CHECKORSAY(size()>0,"assignment with more elements than in array");
    T* first = &(buffer_[0]);
    if (size() > 0)
        *first = e;
    else
        return ra::CommaOp<T>(nullptr, nullptr);
    ra::CommaOp<T> co(first+1, first+size()-1); 
    return co;
}

template<typename T> RA_INLINE_
ra::CommaOp<T>::CommaOp(T* ptr, T* last)
: ptr_(ptr), last_(last)
{ 
    RA_CHECKORSAY(ptr_!=nullptr and last_!=nullptr, "invalid comma operator");
}

template<typename T> RA_INLINE_
ra::CommaOp<T>& ra::CommaOp<T>::operator,(const T& e)
{ 
    RA_CHECKORSAY(ptr_!=nullptr and last_!=nullptr, "invalid comma operator");
    RA_CHECKORSAY(ptr_<=last_, "assignment with more elements than in array");
    if (ptr_ and ptr_ <= last_)
        *ptr_++ = e;
    return *this; 
}


template<typename T, int R> RA_INLINE_
bool ra::rarray<T,R>::is_clear() const
{
    return buffer_.cbegin() == nullptr;
}


template<typename T, int R> RA_INLINE_
ra::rarray<T,R> ra::rarray<T,R>::copy() const {
    rarray<T,R> clone;
    clone.buffer_ = buffer_.copy();
    clone.shape_ = shape_.copy();
    clone.shape_.relocate(clone.buffer_.begin()); // something is wrong with relocate for some types
    return clone;
}


template<typename T, int R> RA_INLINEF
ra::rarray<T,R-1> ra::rarray<T,R>::at(size_type i)
{
    if (i < 0 or i >= extent(0))
        throw std::out_of_range("rarray<T,R>::at");
    size_type stride = size()/extent(0);
    return ra::rarray<T,R-1>(buffer_.slice(i*stride, (i+1)*stride), shape_.at(i));
}

template<typename T, int R> RA_INLINEF
ra::rarray<T,R-1> ra::rarray<T,R>::at(size_type i) const
{
    return const_cast<ra::rarray<T,R>*>(this)->at(i); // provisionary
}

template<typename T, int R> RA_INLINEF
ra::rarray<T,R>::operator typename PointerArray<T,R>::type () 
{
    return shape_.ptrs(); // makes a[..][..] = ... work.
}

template<typename T, int R> RA_INLINEF
ra::rarray<T,R>::operator typename PointerArray<const T,R>::type () const 
{
    return shape_.ptrs(); // makes a[..][..] work.
}

template<typename T, int R> RA_INLINE_
const ra::size_type* ra::rarray<T,R>::shape() const
{
    return &(shape_.extent()[0]);
}

template<typename T, int R> RA_INLINE_
typename ra::rarray<T,R>::parray_t ra::rarray<T,R>::ptr_array() const 
{
    return shape_.ptrs();
}

template<typename T, int R> RA_INLINE_
typename ra::rarray<T,R>::noconst_parray_t ra::rarray<T,R>::noconst_ptr_array() const 
{
    return const_cast<noconst_parray_t>(shape_.ptrs());
}

template<typename T, int R> RA_INLINE_
typename ra::rarray<const T,R>& ra::rarray<T,R>::const_ref() const 
{
    return (rarray<const T,R>&)(*this);
}

template<typename T, int R> RA_INLINEF
ra::rarray<T,R>& ra::rarray<T,R>::operator=(const rarray<T,R> &a) 
{
    buffer_ = a.buffer_;
    shape_ = a.shape_;
    return *this;
}
template<typename T, int R> RA_INLINEF
ra::rarray<T,R>& ra::rarray<T,R>::operator=(rarray<T,R> &&a) 
{
    buffer_ = std::move(a.buffer_);
    shape_ = std::move(a.shape_);
    return *this;
}

template<typename T, int R> RA_INLINEF
void ra::rarray<T,R>::fill(const T& value) {
    buffer_.assign(value);
}

template<typename T, int R> RA_INLINEF
ra::size_type ra::rarray<T,R>::size() const {
    return shape_.size(); 
}

template<typename T, int R> RA_INLINEF
T* ra::rarray<T,R>::data()
{
    return buffer_.begin(); 
}

template<typename T, int R> RA_INLINEF
const T* ra::rarray<T,R>::data() const
{
    return buffer_.begin(); 
}

template<typename T, int R> RA_INLINEF
ra::size_type ra::rarray<T,R>::extent(int i) const
{
    return shape_.extent(i);
}

template<typename T, int R> RA_INLINEF
typename ra::rarray<T,R>::iterator ra::rarray<T,R>::begin()
{
    return buffer_.begin();
}
template<typename T, int R> RA_INLINEF
typename ra::rarray<T,R>::const_iterator ra::rarray<T,R>::begin() const
{
    return buffer_.begin();
}
template<typename T, int R> RA_INLINEF
typename ra::rarray<T,R>::const_iterator ra::rarray<T,R>::cbegin() const
{
    return buffer_.cbegin();
}
template<typename T, int R> RA_INLINEF
typename ra::rarray<T,R>::iterator ra::rarray<T,R>::end()
{
    return buffer_.end();
}
template<typename T, int R> RA_INLINEF
typename ra::rarray<T,R>::const_iterator ra::rarray<T,R>::end() const
{
    return buffer_.end();
}
template<typename T, int R> RA_INLINEF
typename ra::rarray<T,R>::const_iterator ra::rarray<T,R>::cend() const
{
    return buffer_.cend();
}



template<typename T,int R> RA_INLINE_
ra::size_type* ra::rarray<T,R>::index(const iterator&i, size_type* ind) const
{
    return index(*i, ind);
}

template<typename T,int R> RA_INLINE_
ra::size_type* ra::rarray<T,R>::index(const T& a, size_type* ind) const
{
    ptrdiff_t linearindex = &a - &(buffer_[0]);
    RA_CHECKORSAY(linearindex >= 0 and linearindex < size(), "element not in array");
    int j = R;
    const size_type* extent_ = shape();
    while (j-->0) {
        ind[j] = linearindex % extent_[j];
        linearindex /= extent_[j];
    }
    return ind;
}

template<typename T,int R> RA_INLINE_
std::array<ra::size_type,R> ra::rarray<T,R>::index(const iterator&i) const
{
    return index(*i);
}

template<typename T,int R> RA_INLINE_
std::array<ra::size_type,R> ra::rarray<T,R>::index(const T& a) const
{
    std::array<size_type,R> ind;
    ptrdiff_t linearindex = &a - &(buffer_[0]);
    RA_CHECKORSAY(linearindex >=0 and linearindex < size(), "element not in array");
    int j = R;
    const size_type* extent_ = shape();
    while (j-->0) {
        ind[j] = linearindex % extent_[j];
        linearindex /= extent_[j];
    }
    return ind;
}

template<typename T,int R> RA_INLINE_
ra::size_type ra::rarray<T,R>::index(const iterator&iter, int i) const
{
    return index(*iter, i);
}

template<typename T,int R> RA_INLINE_
ra::size_type ra::rarray<T,R>::index(const T& a, int i) const
{
    ptrdiff_t linearindex = &a - &(buffer_[0]);
    if (linearindex < 0 or linearindex >= size())
        throw ("element not in array");
    const size_type* extent_ = shape();
    for (int j = R-1; j > i; j--) 
        linearindex /= extent_[j];
    return linearindex % extent_[i];
}


template<typename A> RA_INLINE_ 
int ra::extent_given_byte_size(A a[], int i, int byte_size) 
{
    RA_CHECKORSAY(i>=0 and i<1, "wrong dimension");
    switch (i) {
        case 0:  return byte_size/sizeof(A);
        default: return 1;
    }
}

template<typename A,int Z> RA_INLINE_ 
int ra::extent_given_byte_size(A a[][Z], int i, int byte_size) 
{
    RA_CHECKORSAY(i>=0 and i<2, "wrong dimension");
    switch (i) {
        case 0:  return byte_size/sizeof(A)/Z;
        case 1:  return Z;
        default: return 1;
    }
}

template<typename A,int Y,int Z> RA_INLINE_ 
int ra::extent_given_byte_size(A a[][Y][Z], int i, int byte_size) 
{
    RA_CHECKORSAY(i>=0 and i<3, "wrong dimension");
    switch (i) {
        case 0:  return byte_size/sizeof(A)/Z/Y;
        case 1:  return Y;
        case 2:  return Z;
        default: return 1;
    }
}

template<typename A,int X,int Y,int Z> RA_INLINE_ 
int ra::extent_given_byte_size(A a[][X][Y][Z], int i, int byte_size) 
{
    RA_CHECKORSAY(i>=0 and i<4, "wrong dimension");
    switch (i) {
        case 0:  return byte_size/sizeof(A)/X/Z/Y;
        case 1:  return X;
        case 2:  return Y;
        case 3:  return Z;
        default: return 1;
    }
}

template<typename A,int W,int X,int Y,int Z> RA_INLINE_ 
int ra::extent_given_byte_size(A a[][W][X][Y][Z], int i, int byte_size) 
{
    RA_CHECKORSAY(i>=0 and i<5, "wrong dimension");
    switch (i) {
        case 0:  return byte_size/sizeof(A)/W/X/Z/Y;
        case 1:  return W;
        case 2:  return X;
        case 3:  return Y;
        case 4:  return Z;
        default: return 1;
    }
}

template<typename A,int V,int W,int X,int Y,int Z> RA_INLINE_ 
int ra::extent_given_byte_size(A a[][V][W][X][Y][Z], int i, int byte_size) 
{
    RA_CHECKORSAY(i>=0 and i<6, "wrong dimension");
    switch (i) {
       case 0:  return byte_size/sizeof(A)/V/W/X/Z/Y;
       case 1:  return V;
       case 2:  return W;
       case 3:  return X;
       case 4:  return Y;
       case 5:  return Z;
       default: return 1;
    }
}

template<typename A,int U,int V,int W,int X,int Y,int Z> RA_INLINE_ 
int ra::extent_given_byte_size(A a[][U][V][W][X][Y][Z], int i, int byte_size) 
{
    RA_CHECKORSAY(i>=0 and i<7, "wrong dimension");
    switch (i) {
       case 0:  return byte_size/sizeof(A)/U/V/W/X/Z/Y;
       case 1:  return U;
       case 2:  return V;
       case 3:  return W;
       case 4:  return X;
       case 5:  return Y;
       case 6:  return Z;
       default: return 1;
    }
}

template<typename A,int T,int U,int V,int W,int X,int Y,int Z> RA_INLINE_ 
int ra::extent_given_byte_size(A a[][T][U][V][W][X][Y][Z], int i, int byte_size) 
{
    RA_CHECKORSAY(i>=0 and i<8, "wrong dimension");
    switch (i) {
       case 0:  return byte_size/sizeof(A)/T/U/V/W/X/Z/Y;
       case 1:  return T;
       case 2:  return U;
       case 3:  return V;
       case 4:  return W;
       case 5:  return X;
       case 6:  return Y;
       case 7:  return Z;
       default: return 1;
    }
}

template<typename A,int S,int T,int U,int V,int W,int X,int Y,int Z> RA_INLINE_ 
int ra::extent_given_byte_size(A a[][S][T][U][V][W][X][Y][Z], int i, int byte_size) 
{
    RA_CHECKORSAY(i>=0 and i<9, "wrong dimension");
    switch (i) {
       case 0:  return byte_size/sizeof(A)/S/T/U/V/W/X/Z/Y;
       case 1:  return S;
       case 2:  return T;
       case 3:  return U;
       case 4:  return V;
       case 5:  return W;
       case 6:  return X;
       case 7:  return Y;
       case 8:  return Z;
       default: return 1;
    }
}

template<typename A,int R,int S,int T,int U,int V,int W,int X,int Y,int Z> RA_INLINE_ 
int ra::extent_given_byte_size(A a[][R][S][T][U][V][W][X][Y][Z], int i, int byte_size) 
{
    RA_CHECKORSAY(i>=0 and i<10, "wrong dimension");
    switch (i) {
       case 0:  return byte_size/sizeof(A)/R/S/T/U/V/W/X/Z/Y;
       case 1:  return R;
       case 2:  return S;
       case 3:  return T;
       case 4:  return U;
       case 5:  return V;
       case 6:  return W;
       case 7:  return X;
       case 8:  return Y;
       case 9:  return Z;
       default: return 1;
    }
}

template<typename A,int Q,int R,int S,int T,int U,int V,int W,int X,int Y,int Z> RA_INLINE_ 
int ra::extent_given_byte_size(A a[][Q][R][S][T][U][V][W][X][Y][Z], int i, int byte_size) 
{
    RA_CHECKORSAY(i>=0 and i<11, "wrong dimension");
    switch (i) {
       case 0:   return byte_size/sizeof(A)/Q/R/S/T/U/V/W/X/Z/Y;
       case 1:   return Q;
       case 2:   return R;
       case 3:   return S;
       case 4:   return T;
       case 5:   return U;
       case 6:   return V;
       case 7:   return W;
       case 8:   return X;
       case 9:   return Y;
       case 10:  return Z;
       default: return 1;
    }
}


template<typename A,int R> RA_INLINE_ 
int ra::extent_given_byte_size(const ra::rarray<A,R>& a, int i, int byte_size) 
{
    return a.extent(i);
}



template<typename A> RA_INLINE_ 
ra::rarray<A,1> ra::make_rarray_given_byte_size(A a[], int byte_size) 
{
    const int z = byte_size/sizeof(A);
    return ra::rarray<A,1>(a,z);
}

template<typename A,int Z> RA_INLINE_ 
ra::rarray<A,2> ra::make_rarray_given_byte_size(A a[][Z], int byte_size) 
{
    const int y = byte_size/sizeof(A)/Z;
    return ra::rarray<A,2>(*a,y,Z);
}

template<typename A,int Y,int Z> RA_INLINE_ 
ra::rarray<A,3> ra::make_rarray_given_byte_size(A a[][Y][Z], int byte_size) 
{
    const int x = byte_size/sizeof(A)/Z/Y;
    return ra::rarray<A,3>(**a,x,Y,Z);
}

template<typename A,int X,int Y,int Z> RA_INLINE_ 
ra::rarray<A,4> ra::make_rarray_given_byte_size(A a[][X][Y][Z], int byte_size) 
{
    const int w = byte_size/sizeof(A)/X/Z/Y;
    return ra::rarray<A,4>(***a,w,X,Y,Z);
}

template<typename A,int W,int X,int Y,int Z> RA_INLINE_ 
ra::rarray<A,5> ra::make_rarray_given_byte_size(A a[][W][X][Y][Z], int byte_size) 
{
    const int v = byte_size/sizeof(A)/W/X/Z/Y;
    return ra::rarray<A,5>(****a,v,W,X,Y,Z);
}

template<typename A,int V,int W,int X,int Y,int Z> RA_INLINE_ 
ra::rarray<A,6> ra::make_rarray_given_byte_size(A a[][V][W][X][Y][Z], int byte_size) 
{
    const int u = byte_size/sizeof(A)/V/W/X/Z/Y;
    return ra::rarray<A,6>(*****a,u,V,W,X,Y,Z);
}

template<typename A,int U,int V,int W,int X,int Y,int Z> RA_INLINE_ 
ra::rarray<A,7> ra::make_rarray_given_byte_size(A a[][U][V][W][X][Y][Z], int byte_size) 
{
    const int t = byte_size/sizeof(A)/U/V/W/X/Z/Y;
    return ra::rarray<A,7>(******a,t,U,V,W,X,Y,Z);
}

template<typename A,int T,int U,int V,int W,int X,int Y,int Z> RA_INLINE_ 
ra::rarray<A,8> ra::make_rarray_given_byte_size(A a[][T][U][V][W][X][Y][Z], int byte_size) 
{
    const int s = byte_size/sizeof(A)/T/U/V/W/X/Z/Y;
    return ra::rarray<A,8>(*******a,s,T,U,V,W,X,Y,Z);
}

template<typename A,int S,int T,int U,int V,int W,int X,int Y,int Z> RA_INLINE_ 
ra::rarray<A,9> ra::make_rarray_given_byte_size(A a[][S][T][U][V][W][X][Y][Z], int byte_size) 
{
    const int r = byte_size/sizeof(A)/S/T/U/V/W/X/Z/Y;
    return ra::rarray<A,9>(********a,r,S,T,U,V,W,X,Y,Z);
}

template<typename A,int R,int S,int T,int U,int V,int W,int X,int Y,int Z> RA_INLINE_ 
ra::rarray<A,10> ra::make_rarray_given_byte_size(A a[][R][S][T][U][V][W][X][Y][Z], int byte_size) 
{
    const int q = byte_size/sizeof(A)/R/S/T/U/V/W/X/Z/Y;
    return ra::rarray<A,10>(*********a,q,R,S,T,U,V,W,X,Y,Z);
}

template<typename A,int Q,int R,int S,int T,int U,int V,int W,int X,int Y,int Z> RA_INLINE_ 
ra::rarray<A,11> ra::make_rarray_given_byte_size(A a[][Q][R][S][T][U][V][W][X][Y][Z], int byte_size) 
{
    const int p = byte_size/sizeof(A)/Q/R/S/T/U/V/W/X/Z/Y;
    return ra::rarray<A,11>(**********a,p,Q,R,S,T,U,V,W,X,Y,Z);
}


template<typename T,int R> RA_INLINE_ 
ra::rarray<T,R> ra::make_rarray_given_byte_size(ra::rarray<T,R> a, int byte_size) 
{
    return a;
}



#undef RA_CHECKORSAY
#undef RA_INLINE
#undef RA_INLINEF
#undef RA_INLINE_



    
#define EXTENT(A,I)  ra::extent_given_byte_size(A,I,sizeof(A))
#define RARRAY(A)    ra::make_rarray_given_byte_size(A,sizeof(A))
#define INDEX(A,X,I) RARRAY(A).index(X,I)
#define as_rarray(A) ra::make_rarray_given_byte_size(A,sizeof(A)) 

using ra::rarray;
using ra::linspace;
using ra::xrange; 

template<typename T> using rvector = rarray<T,1>;
template<typename T> using rmatrix = rarray<T,2>;
template<typename T> using rtensor = rarray<T,3>;

#else
#error "This file requires C++11 or newer support."
#endif

#endif