//
// rarray - Runtime arrays: template classes for pointer-based,
//          runtime, reference counted, multi-dimensional
//          arrays.  Documentation in rarraydoc.pdf
//
// Copyright (c) 2013-2022  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

//begin #include "versionheader.h"
#define RA_VERSION "v2.5.0"
#define RA_VERSION_NUMBER 2005000

//end #include "versionheader.h"

//begin #include "rarraymacros.h"
#ifndef _RARRAYMACROSH_
#define _RARRAYMACROSH_

#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__) + " (rarray:" + std::to_string(__LINE__) + ")")
#define RA_NOEXCEPT(x)
#define noboundscheck false
#else
#define RA_CHECKORSAY(a, b)
#define RA_NOEXCEPT(x) noexcept(x)
#define noboundscheck true
#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

#endif

//end #include "rarraymacros.h"

//begin #include "shared_buffer.h"


#ifndef SHAREDBUFFERH
#define SHAREDBUFFERH

//#include "rarraymacros.h" was already done above
#include <cstddef>
#include <algorithm>
#include <stdexcept>
#include <memory>
#include <cstring>
#include <atomic>

namespace ra { template<class T> class shared_buffer; }
template<typename V>
int internal_check(const ra::shared_buffer<V>& a,
                   bool datavalue_shouldbe, V* datavalue,
                   bool origvalue_shouldbe, V* origvalue,
                   bool refsvalue_shouldbe, std::atomic<int>* refsvalue,
                   bool refscount_shouldbe, int refscount,
                   bool sizevalue_shouldbe, typename ra::shared_buffer<V>::size_type sizevalue);

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

namespace ra {

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

template<class T>
class shared_buffer
{
  public:

    typedef ssize_t size_type;

    shared_buffer() noexcept;
    explicit shared_buffer(size_type asize);
    shared_buffer(size_type asize, T* adata) RA_NOEXCEPT(true);
    shared_buffer(const shared_buffer& other) noexcept;
    shared_buffer(shared_buffer&& from) noexcept;

    shared_buffer& operator=(const shared_buffer& other) noexcept;
    void operator=(shared_buffer&& from) noexcept;

    ~shared_buffer() noexcept;

    const T& operator[](size_type index) const RA_NOEXCEPT(true);
    T& operator[](size_type index) RA_NOEXCEPT(true);

    const T& at(size_type index) const;
    T& at(size_type index);

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

    size_type size() const noexcept;

    shared_buffer<T> copy() const;

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

    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_;

    template<typename V>
    friend int ::internal_check(const ra::shared_buffer<V>& a,
                   bool datavalue_shouldbe, V* datavalue,
                   bool origvalue_shouldbe, V* origvalue,
                   bool refsvalue_shouldbe, std::atomic<int>* refsvalue,
                   bool refscount_shouldbe, int refscount,
                   bool sizevalue_shouldbe, typename ra::shared_buffer<V>::size_type sizevalue);
    void uninit() noexcept;
    void incref() noexcept;
    void decref() noexcept;
};

/***************************************************************************/
/* FUTURE FEATURE :
#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 {
            free(some_orig);
            throw std::bad_alloc();
        }
    } else {
        throw std::bad_alloc();
    }
}
*/

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

template<class T>
shared_buffer<T>::shared_buffer(size_type asize)
  : data_(nullptr), orig_(nullptr), size_(0), refs_(nullptr)
{
    T* data;
    refs_ = new std::atomic<int>(1);
    try {
        data = new T[asize];
    }
    catch (...) {
        delete refs_;
        throw;
    }
    data_ = data;
    orig_ = data_;
    size_ = asize;
}

template<class T>
shared_buffer<T>::shared_buffer(size_type asize, T* adata) RA_NOEXCEPT(true)
  : data_(adata), orig_(nullptr), size_(asize), refs_(nullptr)
{
    RA_CHECKORSAY(adata, "nullptr given as data");
}

template<class T>
shared_buffer<T>::shared_buffer(const shared_buffer<T>& other) noexcept
  : 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) noexcept
  : 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) noexcept
{
    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) noexcept
{
    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 RA_NOEXCEPT(true)
{
    RA_CHECKORSAY(index >= 0 and index < size(), "element not in buffer");
    return data_[index];
}

template<class T>
T& shared_buffer<T>::operator[](size_type index) RA_NOEXCEPT(true)
{
    RA_CHECKORSAY(index >= 0 and index < size(), "element not in buffer");
    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>
const shared_buffer<T> shared_buffer<T>::slice(size_type from, size_type to) const
{
    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() noexcept
{
    decref();
}

template<class T>
typename shared_buffer<T>::size_type shared_buffer<T>::size() const noexcept
{
    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() noexcept
{
    return data_;
}

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

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

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

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

template<class T>
typename shared_buffer<T>::const_iterator shared_buffer<T>::end() const noexcept
{
    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;
    } else {
        std::atomic<int>* newrefs = new std::atomic<int>(1);
        T* newdata;
        try {
             newdata = new T[newsize];
        }
        catch (...) {
            delete newrefs;
            throw;
        }
        if (keep_content) {
            try {
                size_t n = ((size_<newsize)?size_:newsize);
                for (size_t i=0; i<n; i++)
                    newdata[i] = data_[i];
            }
            catch (...) {
                delete newrefs;
                delete[] newdata;
                throw;
            }
        }
        decref();
        data_ = newdata;
        orig_ = newdata;
        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() noexcept {
    data_ = nullptr;
    orig_ = nullptr;
    size_ = 0;
    refs_ = nullptr;
}

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

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

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

//end #include "shared_buffer.h"

//begin #include "shared_shape.h"

#ifndef SHARED_SHAPEH
#define SHARED_SHAPEH

//#include "rarraymacros.h" was already done above
//begin #include "offsets.h"
#ifndef OFFSETSH
#define OFFSETSH

//#include "rarraymacros.h" was already done above
#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 noexcept;
    ssize_t get_num_offsets() const noexcept;
    ssize_t get_rank() const noexcept;

  private:

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

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

inline Offsets::Offsets(const std::vector<ssize_t>& extent)
{
    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) {
            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
{
    static_assert(sizeof(T*) == sizeof(void*) &&
                  sizeof(T*) == sizeof(void**),
                  "rarray's Offsets requires all pointers to have the same size");
    ssize_t noffsets = offsets_.size();
    if (ndataoffsets_ == 0 && noffsets == 0)
        return nullptr;
    else if (ndataoffsets_ == 1 && noffsets == 0)
        return reinterpret_cast<void***>(data);
    else {
        void*** offsets = new void**[noffsets];
        ssize_t i = 0;
        for (;i < noffsets - ndataoffsets_; i++)
            offsets[i] = reinterpret_cast<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 noexcept
{
    return ndataoffsets_;
}

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

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

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

#endif

//end #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;
    typedef typename PointerArray<T,R-1>::noconst_type*  noconst_type;
};
template<typename T>
struct PointerArray<T,1> {
    typedef T* type;
    typedef T* noconst_type;
};
template<typename T>
struct PointerArray<T,0> {
    typedef T& type;
    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() noexcept;
    shared_shape(const std::array<size_type,R>&anextent, T*adata);
    shared_shape(const shared_shape& other) noexcept;
    shared_shape(shared_shape&& other) noexcept;

    shared_shape& operator=(const shared_shape& other) noexcept;
    void operator=(shared_shape&& other) noexcept;

    ~shared_shape();

    shared_shape copy() const;

    void relocate(T* newdata);
    void reshape(const std::array<size_type,R>&newextent);

    ptrs_type ptrs() const noexcept;
    T* data() const noexcept;
    size_type size() const noexcept;
    size_type extent(int i) const;
    const std::array<size_type,R>& extent() const noexcept;

    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() noexcept;
    void incref() const noexcept;
    void decref() noexcept;
    void copy_before_write();

    template<class U, int S> friend class shared_shape;

    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() const noexcept {}
    template<class U, int S> friend class shared_shape;
};

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

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

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

template<class T, int R>
shared_shape<T,R>::shared_shape(const shared_shape& other) noexcept
  : 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) noexcept
  : 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) noexcept
{
    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) noexcept
{
    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) noexcept
    {
        return reinterpret_cast<T*>(
                  const_cast<typename PointerArray<T,R-1>::noconst_type>(
                    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) noexcept
    {
        return reinterpret_cast<T*>(ptrs);
    }
};

template<class T, int R>
T* shared_shape<T,R>::data() const noexcept
{
    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 noexcept
{
    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
{
    if (i < 0 or i >= R)
        throw std::out_of_range("shared_shape::extent(int)");
    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 noexcept {
    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 = reinterpret_cast<char*>(newdata) - reinterpret_cast<char*>(data());
        if (shift != 0) {
            copy_before_write();
            for (size_type i = noffsets_ - ndataoffsets_; i < noffsets_; i++)
                orig_[i] = reinterpret_cast<void**>(reinterpret_cast<char*>(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_];
        void*** old_eff_orig = reinterpret_cast<void***>(
                                   const_cast<typename PointerArray<T,R>::noconst_type>(
                                     ptrs_ ));
        std::copy(old_eff_orig, old_eff_orig + noffsets_, copy_of_this.orig_);
        std::ptrdiff_t shift = reinterpret_cast<char*>(copy_of_this.orig_)
                               - reinterpret_cast<char*>(old_eff_orig);
        for (size_type i = 0; i < noffsets_ - ndataoffsets_; i++)
            copy_of_this.orig_[i] = reinterpret_cast<void**>(reinterpret_cast<char*>(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>&newextent)
{
    if (newextent != extent_) {
        *this = shared_shape<T,R>(newextent, data());
    }
}

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

template<class T, int R>
void shared_shape<T,R>::uninit() noexcept
{
    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 noexcept
{
    if (refs_)
        (*refs_)++;
}

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

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

}

#endif

//end #include "shared_shape.h"


#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <list>
#include <sstream>
#include <stdexcept>
#include <string>
#include <utility>

namespace ra {
typedef ssize_t size_type;
enum class RESIZE { NO, ALLOWED };
template<typename T,int R> class rarray;
template<typename T> class CommaOp;
template<typename T,int R> RA_INLINE_ std::istream& operator>>(std::istream &i, rarray<T,R>& r);
template<typename T,int R> RA_INLINE_ std::ostream& operator<<(std::ostream &o, const rarray<T,R>& r);
template<typename T,int R,typename P> class _Bracket;
template<typename T,int R,typename P> class _ConstBracket;
#define ExOp class
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 T                                         value_type;
    typedef ssize_t                                   difference_type;
    typedef ssize_t                                   size_type;
    typedef T*                                        iterator;
    typedef const T*                                  const_iterator;
    typedef typename PointerArray<T,R>::type          parray_t;
    typedef typename PointerArray<T,R>::noconst_type  noconst_parray_t;
    RA_INLINE_ rarray();
    template<int R_=R,class=typename std::enable_if<R_==1>::type>
    RA_INLINE_ explicit rarray(size_type n0);
    template<int R_=R,class=typename std::enable_if<R_==2>::type>
    RA_INLINE_ rarray(size_type n0, size_type n1);
    template<int R_=R,class=typename std::enable_if<R_==3>::type>
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2);
    template<int R_=R,class=typename std::enable_if<R_==4>::type>
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2, size_type n3);
    template<int R_=R,class=typename std::enable_if<R_==5>::type>
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4);
    template<int R_=R,class=typename std::enable_if<R_==6>::type>
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5);
    template<int R_=R,class=typename std::enable_if<R_==7>::type>
    RA_INLINE_ rarray(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5, size_type n6);
    template<int R_=R,class=typename std::enable_if<R_==8>::type>
    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);
    template<int R_=R,class=typename std::enable_if<R_==9>::type>
    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);
    template<int R_=R,class=typename std::enable_if<R_==10>::type>
    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);
    template<int R_=R,class=typename std::enable_if<R_==11>::type>
    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);
    RA_INLINE_ rarray(const size_type* anextent);
    template<int R_=R,class=typename std::enable_if<R_==1>::type>
    RA_INLINE_ rarray(T* buffer, size_type n0);
    template<int R_=R,class=typename std::enable_if<R_==2>::type>
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1);
    template<int R_=R,class=typename std::enable_if<R_==3>::type>
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1, size_type n2);
    template<int R_=R,class=typename std::enable_if<R_==4>::type>
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1, size_type n2, size_type n3);
    template<int R_=R,class=typename std::enable_if<R_==5>::type>
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1, size_type n2, size_type n3, size_type n4);
    template<int R_=R,class=typename std::enable_if<R_==6>::type>
    RA_INLINE_ rarray(T* buffer, size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, size_type n5);
    template<int R_=R,class=typename std::enable_if<R_==7>::type>
    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);
    template<int R_=R,class=typename std::enable_if<R_==8>::type>
    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);
    template<int R_=R,class=typename std::enable_if<R_==9>::type>
    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);
    template<int R_=R,class=typename std::enable_if<R_==10>::type>
    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);
    template<int R_=R,class=typename std::enable_if<R_==11>::type>
    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);
    RA_INLINE_ rarray(T* buffer, const size_type* anextent);
    template<std::size_t Z,
             int R__=R,typename=typename std::enable_if<R__==1>::type>
    RA_INLINE_ explicit rarray(T (&a)[Z]);
    template<std::size_t Y,std::size_t Z,
             int R__=R,typename=typename std::enable_if<R__==2>::type>
    RA_INLINE_ explicit rarray(T (&a)[Y][Z]);
    template<std::size_t X,std::size_t Y,std::size_t Z,
             int R__=R,typename=typename std::enable_if<R__==3>::type>
    RA_INLINE_ explicit rarray(T (&a)[X][Y][Z]);
    template<std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,
             int R__=R,typename=typename std::enable_if<R__==4>::type>
    RA_INLINE_ explicit rarray(T (&a)[W][X][Y][Z]);
    template<std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,
             int R__=R,typename=typename std::enable_if<R__==5>::type>
    RA_INLINE_ explicit rarray(T (&a)[V][W][X][Y][Z]);
    template<std::size_t U,std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,
             int R__=R,typename=typename std::enable_if<R__==6>::type>
    RA_INLINE_ explicit rarray(T (&a)[U][V][W][X][Y][Z]);
    template<std::size_t T_,std::size_t U,std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,
             int R__=R,typename=typename std::enable_if<R__==7>::type>
    RA_INLINE_ explicit rarray(T (&a)[T_][U][V][W][X][Y][Z]);
    template<std::size_t S,std::size_t T_,std::size_t U,std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,
             int R__=R,typename=typename std::enable_if<R__==8>::type>
    RA_INLINE_ explicit rarray(T (&a)[S][T_][U][V][W][X][Y][Z]);
    template<std::size_t R_,std::size_t S,std::size_t T_,std::size_t U,std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,
             int R__=R,typename=typename std::enable_if<R__==9>::type>
    RA_INLINE_ explicit rarray(T (&a)[R_][S][T_][U][V][W][X][Y][Z]);
    template<std::size_t Q,std::size_t R_,std::size_t S,std::size_t T_,std::size_t U,std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,
             int R__=R,typename=typename std::enable_if<R__==10>::type>
    RA_INLINE_ explicit rarray(T (&a)[Q][R_][S][T_][U][V][W][X][Y][Z]);
    template<std::size_t P,std::size_t Q,std::size_t R_,std::size_t S,std::size_t T_,std::size_t U,std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,
             int R__=R,typename=typename std::enable_if<R__==11>::type>
    RA_INLINE_ explicit rarray(T (&a)[P][Q][R_][S][T_][U][V][W][X][Y][Z]);
    RA_INLINEF rarray(const rarray<T,R> &a) noexcept;
    RA_INLINE_ rarray<T,R>& operator=(const rarray<T,R> &a) noexcept;
    rarray(rarray<T,R>&& x) noexcept;
    rarray<T,R>& operator=(rarray<T,R>&& x) noexcept;
    RA_INLINE_ CommaOp<T> operator=(const T& e) RA_NOEXCEPT(std::is_nothrow_copy_constructible<T>());
    RA_INLINEF ~rarray();
    constexpr int rank() const noexcept { return R; }
    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() noexcept;
    template<int R_=R,class=typename std::enable_if<R_==1>::type>
    RA_INLINE_ void reshape(size_type n0, RESIZE resize_allowed=RESIZE::NO);
    template<int R_=R,class=typename std::enable_if<R_==2>::type>
    RA_INLINE_ void reshape(size_type n0, size_type n1, RESIZE resize_allowed=RESIZE::NO);
    template<int R_=R,class=typename std::enable_if<R_==3>::type>
    RA_INLINE_ void reshape(size_type n0, size_type n1, size_type n2, RESIZE resize_allowed=RESIZE::NO);
    template<int R_=R,class=typename std::enable_if<R_==4>::type>
    RA_INLINE_ void reshape(size_type n0, size_type n1, size_type n2, size_type n3, RESIZE resize_allowed=RESIZE::NO);
    template<int R_=R,class=typename std::enable_if<R_==5>::type>
    RA_INLINE_ void reshape(size_type n0, size_type n1, size_type n2, size_type n3, size_type n4, RESIZE resize_allowed=RESIZE::NO);
    template<int R_=R,class=typename std::enable_if<R_==6>::type>
    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);
    template<int R_=R,class=typename std::enable_if<R_==7>::type>
    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);
    template<int R_=R,class=typename std::enable_if<R_==8>::type>
    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);
    template<int R_=R,class=typename std::enable_if<R_==9>::type>
    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);
    template<int R_=R,class=typename std::enable_if<R_==10>::type>
    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);
    template<int R_=R,class=typename std::enable_if<R_==11>::type>
    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);
    RA_INLINE_ void reshape(const size_type* extent, RESIZE resize_allowed=RESIZE::NO);
    RA_INLINE_ bool                empty()           const noexcept;
    RA_INLINE_ rarray<T,R>         copy()               const;
    RA_INLINE_ size_type           extent(int i)        const;
    RA_INLINEF const size_type*    shape()              const noexcept;
    RA_INLINE_ size_type           size()               const noexcept;
    RA_INLINE_ T*                  data()               noexcept;
    RA_INLINE_ const T*            data()               const noexcept;
    RA_INLINE_ parray_t            ptr_array()          const noexcept;
    RA_INLINE_ noconst_parray_t    noconst_ptr_array()  const noexcept;
    RA_INLINE_ const rarray<const T,R>&  const_ref()    const noexcept;
    RA_INLINE_ void                fill(const T& value);
    RA_INLINE_ iterator            begin()              noexcept;
    RA_INLINE_ const_iterator      begin()              const noexcept;
    RA_INLINE_ const_iterator      cbegin()             const noexcept;
    RA_INLINE_ iterator            end()                noexcept;
    RA_INLINE_ const_iterator      end()                const noexcept;
    RA_INLINE_ const_iterator      cend()               const noexcept;
    RA_INLINE_ size_type           index(const T& a, int i) const;
    RA_INLINE_ size_type           index(const iterator& iter, int i) const;
    RA_INLINE_ std::array<size_type,R> index(const T& a) const;
    RA_INLINE_ std::array<size_type,R> index(const iterator& i) const;
    RA_INLINE_ size_type*          index(const T& a, size_type* ind) const;
    RA_INLINE_ size_type*          index(const iterator& i, size_type* ind) const;
    template<class U=rarray<T,R-1>>
    RA_INLINEF typename std::enable_if<R!=1,U>::type at(size_type i);
    template<class U=T&>
    RA_INLINEF typename std::enable_if<R==1,U>::type at(size_type i);
    template<class U=const rarray<T,R-1>>
    RA_INLINEF typename std::enable_if<R!=1,U>::type at(size_type i) const;
    template<class U=const T&>
    RA_INLINEF typename std::enable_if<R==1,U>::type at(size_type i) const;
    /* * *
    RA_INLINEF operator typename PointerArray<T,R>::type () noexcept;
    RA_INLINEF operator typename PointerArray<const T,R>::type () const noexcept ;
    * * */
    template<class U=T> typename std::enable_if<R==1,
    U&
    >::type
    operator[](ssize_t index) {
        RA_CHECKORSAY(index >= 0 and index < extent(0), "index out of range of array");
        return buffer_.begin()[index];
    }
    template<class U=T> typename std::enable_if<R!=1,
    typename ra::_Bracket<U,R-1,rarray>
    >::type
    operator[](ssize_t index) {
        return { *this, index, this->shape() };
    }
    template<class U=T> typename std::enable_if<R==1,
    const U&
    >::type
    operator[](ssize_t index) const {
        RA_CHECKORSAY(index >= 0 and index < extent(0), "index out of range of array");
        return buffer_.cbegin()[index];
    }
    template<class U=T> typename std::enable_if<R!=1,
    typename ra::_ConstBracket<U,R-1,rarray>
    >::type
    operator[](ssize_t index) const {
        return { *this, index, this->shape() };
    }
    RA_INLINEF const T& leval(size_type i) const;
  private:
    friend class _Bracket<T,R-1,rarray<T,R>>;
    friend class _ConstBracket<T,R-1,rarray<T,R>>;
    RA_INLINEF typename PointerArray<T,R-1>::type _at(ssize_t index);
    RA_INLINEF typename PointerArray<T,R-1>::type _at(ssize_t index) const;
    friend class rarray<T,R+1>;
    rarray(shared_buffer<T>&& abuffer, shared_shape<T,R>&& ashape) :
        buffer_(std::forward<shared_buffer<T>>(abuffer)),
        shape_(std::forward<shared_shape<T,R>>(ashape))
    {}
    rarray(const shared_buffer<T>&& abuffer, shared_shape<T,R>&& ashape) :
        buffer_(std::forward<shared_buffer<T>>(const_cast<shared_buffer<T>&& >(abuffer))),
        shape_(std::forward<shared_shape<T,R>>(ashape))
    {}
    shared_buffer<T>  buffer_;
    shared_shape<T,R> shape_;

};

template<typename T>
class CommaOp {
  public:
    RA_INLINEF CommaOp& operator,(const T& e);
  private:
    RA_INLINEF CommaOp(T* ptr, T* last) RA_NOEXCEPT(true);
    T *ptr_;
    T * const last_;
    template<typename,int> friend class rarray;
    template<typename,int> friend class subrarray;
};

template<typename T,int R,typename P>
class _Bracket {
  private:
    P&             parent_; // what array/array expre is being accessed
    ssize_t        index_;
    const ssize_t* shape_;
  public:
    RA_INLINEF _Bracket<T,R-1,_Bracket> operator[](ssize_t nextindex) noexcept(noboundscheck);
    RA_INLINEF operator decltype(parent_.at(index_)) ();
    _Bracket(const _Bracket&) = delete;
    _Bracket(const _Bracket&&) = delete;
    _Bracket& operator=(const _Bracket&) = delete;
    _Bracket& operator=(const _Bracket&&) = delete;
  private:
    RA_INLINEF _Bracket(P& parent, ssize_t index, const ssize_t* shape) noexcept(noboundscheck);
    RA_INLINEF auto _at(ssize_t nextindex) noexcept(noboundscheck) -> decltype(parent_._at(index_)[nextindex]);
    template<typename T2,int R2> friend class ra::rarray;
    template<typename T2,int R2,typename P2> friend class ra::_Bracket;
  public:
    RA_INLINEF auto at(ssize_t nextindex) -> decltype(parent_.at(index_).at(nextindex));
};

template<typename T,typename P>
class _Bracket<T,1,P> {
  private:
    P&             parent_; // what array/array expre is being accessed
    ssize_t        index_;
    const ssize_t* shape_;
  public:
    RA_INLINEF T& operator[](ssize_t nextindex) noexcept(noboundscheck);
    RA_INLINEF operator decltype(parent_.at(index_)) ();
    _Bracket(const _Bracket&) = delete;
    _Bracket(const _Bracket&&) = delete;
    _Bracket& operator=(const _Bracket&) = delete;
    _Bracket& operator=(const _Bracket&&) = delete;
  private:
    RA_INLINEF _Bracket(P& parent, ssize_t index, const ssize_t* shape) noexcept(noboundscheck);
    RA_INLINEF auto _at(ssize_t nextindex) noexcept(noboundscheck) -> decltype(parent_._at(index_)[nextindex]);
    template<typename T2,int R2> friend class ra::rarray;
    template<typename T2,int R2,typename P2> friend class ra::_Bracket;
  public:
    RA_INLINEF auto at(ssize_t nextindex) -> decltype(parent_.at(index_).at(nextindex));
};

template<typename T,int R,typename P>
class _ConstBracket {
  private:
    const P&       parent_;
    ssize_t        index_;
  public:
    RA_INLINEF _ConstBracket<T,R-1,_ConstBracket> operator[](ssize_t nextindex) const noexcept(noboundscheck);
    RA_INLINEF operator decltype(parent_.at(index_)) ();
    _ConstBracket(const _ConstBracket&) = delete;
    _ConstBracket(const _ConstBracket&&) = delete;
    _ConstBracket& operator=(const _ConstBracket&) = delete;
    _ConstBracket& operator=(const _ConstBracket&&) = delete;
  private:
    const ssize_t* shape_;
    RA_INLINEF _ConstBracket(const P& parent, ssize_t index, const ssize_t* shape) noexcept(noboundscheck);
    RA_INLINEF auto _at(ssize_t nextindex) const noexcept(noboundscheck) -> decltype(parent_._at(index_)[nextindex]);
    template<typename T2,int R2> friend class ra::rarray;
    template<typename T2,int R2,typename P2> friend class ra::_ConstBracket;
  public:
    RA_INLINEF auto at(ssize_t nextindex) const -> decltype(parent_.at(index_).at(nextindex));
};

template<typename T,typename P>
class _ConstBracket<T,1,P> {
  private:
    const P&       parent_;
    ssize_t        index_;
  public:
    RA_INLINEF const T& operator[](ssize_t nextindex) const noexcept(noboundscheck);
    RA_INLINEF operator decltype(parent_.at(index_)) ();
    _ConstBracket(const _ConstBracket&) = delete;
    _ConstBracket(const _ConstBracket&&) = delete;
    _ConstBracket& operator=(const _ConstBracket&) = delete;
    _ConstBracket& operator=(const _ConstBracket&&) = delete;
  private:
    const ssize_t* shape_;
    RA_INLINEF _ConstBracket(const P& parent, ssize_t index, const ssize_t* shape) noexcept(noboundscheck);
    RA_INLINEF auto _at(ssize_t nextindex) const noexcept(noboundscheck) -> decltype(parent_._at(index_)[nextindex]);
    template<typename T2,int R2> friend class ra::rarray;
    template<typename T2,int R2,typename P2> friend class ra::_ConstBracket;
  public:
    RA_INLINEF auto at(ssize_t nextindex) const -> decltype(parent_.at(index_).at(nextindex));
};

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 = static_cast<int>(x2 - x1 + end_incl);
        else
            n = static_cast<int>(x1 - x2 + end_incl);
    }
    rarray<S,1> x(n);
    for (int i = 0; i < n; i++)
        x[i] = x1 + static_cast<S>(((x2-x1)*static_cast<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 + static_cast<T>(static_cast<T>(std::ceil(static_cast<double>(b-a)/static_cast<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>(static_cast<T>(0), end, static_cast<T>(1));
}

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

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

}


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

template<typename T, int R>
template<int R_,typename>
RA_INLINE_
ra::rarray<T,R>::rarray(ra::size_type n0)
  : buffer_(n0),
    shape_({n0}, buffer_.begin())
{}

template<typename T, int R>
template<int R_,typename>
RA_INLINE_
ra::rarray<T,R>::rarray(size_type n0,size_type n1)
  : buffer_(n0*n1),
    shape_({n0,n1}, buffer_.begin())
{}

template<typename T, int R>
template<int R_,typename>
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())
{}

template<typename T, int R>
template<int R_,typename>
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())
{}

template<typename T, int R>
template<int R_,typename>
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())
{}

template<typename T, int R>
template<int R_,typename>
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())
{}

template<typename T, int R>
template<int R_,typename>
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())
{}

template<typename T, int R>
template<int R_,typename>
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())
{}

template<typename T, int R>
template<int R_,typename>
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())
{}

template<typename T, int R>
template<int R_,typename>
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())
{}

template<typename T, int R>
template<int R_,typename>
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())
{}

template<typename T, int R>
template<std::size_t Z,int R__,typename> RA_INLINE_
ra::rarray<T,R>::rarray(T (&a)[Z])
  : buffer_(Z,  a),
    shape_({Z}, buffer_.begin())
{}

template<typename T, int R>
template<std::size_t Y,std::size_t Z,int R__,typename> RA_INLINE_
ra::rarray<T,R>::rarray(T (&a)[Y][Z])
   : buffer_(Y*Z,  *a),
     shape_({Y,Z}, buffer_.begin())
{}

template<typename T, int R>
template<std::size_t X,std::size_t Y,std::size_t Z,int R__,typename> RA_INLINE_
ra::rarray<T,R>::rarray(T (&a)[X][Y][Z])
   : buffer_(X*Y*Z,  **a),
     shape_({X,Y,Z}, buffer_.begin())
{}

template<typename T, int R>
template<std::size_t W, std::size_t X,std::size_t Y,std::size_t Z,int R__,typename> RA_INLINE_
ra::rarray<T,R>::rarray(T (&a)[W][X][Y][Z])
   : buffer_(W*X*Y*Z,  ***a),
     shape_({W,X,Y,Z}, buffer_.begin())
{}

template<typename T, int R>
template<std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,int R__,typename> RA_INLINE_
ra::rarray<T,R>::rarray(T (&a)[V][W][X][Y][Z])
   : buffer_(V*W*X*Y*Z,  ****a),
     shape_({V,W,X,Y,Z}, buffer_.begin())
{}

template<typename T, int R>
template<std::size_t U,std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,int R__,typename> RA_INLINE_
ra::rarray<T,R>::rarray(T (&a)[U][V][W][X][Y][Z])
   : buffer_(U*V*W*X*Y*Z,  *****a),
     shape_({U,V,W,X,Y,Z}, buffer_.begin())
{}

template<typename T, int R>
template<std::size_t T_,std::size_t U,std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,int R__,typename> RA_INLINE_
ra::rarray<T,R>::rarray(T (&a)[T_][U][V][W][X][Y][Z])
   : buffer_(T_*U*V*W*X*Y*Z,  ******a),
     shape_({T_,U,V,W,X,Y,Z}, buffer_.begin())
{}

template<typename T, int R>
template<std::size_t S,std::size_t T_,std::size_t U,std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,int R__,typename> RA_INLINE_
ra::rarray<T,R>::rarray(T (&a)[S][T_][U][V][W][X][Y][Z])
   : buffer_(S*T_*U*V*W*X*Y*Z,  *******a),
     shape_({S,T_,U,V,W,X,Y,Z}, buffer_.begin())
{}

template<typename T, int R>
template<std::size_t R_,std::size_t S,std::size_t T_,std::size_t U,std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,int R__,typename> RA_INLINE_
ra::rarray<T,R>::rarray(T (&a)[R_][S][T_][U][V][W][X][Y][Z])
   : buffer_(R_*S*T_*U*V*W*X*Y*Z,  ********a),
     shape_({R_,S,T_,U,V,W,X,Y,Z}, buffer_.begin())
{}

template<typename T, int R>
template<std::size_t Q,std::size_t R_,std::size_t S,std::size_t T_,std::size_t U,std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,int R__,typename> RA_INLINE_
ra::rarray<T,R>::rarray(T (&a)[Q][R_][S][T_][U][V][W][X][Y][Z])
   : buffer_(Q*R_*S*T_*U*V*W*X*Y*Z,  *********a),
     shape_({Q,R_,S,T_,U,V,W,X,Y,Z}, buffer_.begin())
{}

template<typename T, int R>
template<std::size_t P,std::size_t Q,std::size_t R_,std::size_t S,std::size_t T_,std::size_t U,std::size_t V,std::size_t W,std::size_t X,std::size_t Y,std::size_t Z,int R__,typename> RA_INLINE_
ra::rarray<T,R>::rarray(T (&a)[P][Q][R_][S][T_][U][V][W][X][Y][Z])
   : buffer_(P*Q*R_*S*T_*U*V*W*X*Y*Z,  **********a),
     shape_({P,Q,R_,S,T_,U,V,W,X,Y,Z}, buffer_.begin())
{}

template<typename T, int R>
template<int R_,typename>
RA_INLINE_
ra::rarray<T,R>::rarray(T* buffer, size_type n0)
  : buffer_(n0,  buffer),
    shape_({n0}, buffer)
{}

template<typename T, int R>
template<int R_,typename>
RA_INLINE_
ra::rarray<T,R>::rarray(T* buffer, size_type n0,size_type n1)
  : buffer_(n0*n1,  buffer),
    shape_({n0,n1}, buffer)
{}

template<typename T, int R>
template<int R_,typename>
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)
{}

template<typename T, int R>
template<int R_,typename>
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)
{}

template<typename T, int R>
template<int R_,typename>
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)
{}

template<typename T, int R>
template<int R_,typename>
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)
{}

template<typename T, int R>
template<int R_,typename>
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)
{}

template<typename T, int R>
template<int R_,typename>
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)
{}

template<typename T, int R>
template<int R_,typename>
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)
{}

template<typename T, int R>
template<int R_,typename>
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)
{}

template<typename T, int R>
template<int R_,typename>
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)
{}

namespace ra {
RA_INLINE_ size_type mul(const size_type * x, std::size_t n) noexcept {
    size_type result = 1;
    for (std::size_t 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* anextent)
  : buffer_(mul(anextent,R)),
    shape_(reinterpret_cast<const std::array<size_type,R>&>(*anextent), buffer_.begin())
{}

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

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

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


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


template<typename T, int R>
template<int R_,typename>
RA_INLINE_
void ra::rarray<T,R>::reshape(size_type n0, RESIZE resize_allowed)
{
    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>
template<int R_,typename>
RA_INLINE_
void ra::rarray<T,R>::reshape(size_type n0, size_type n1, RESIZE resize_allowed)
{
    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>
template<int R_,typename>
RA_INLINE_
void ra::rarray<T,R>::reshape(size_type n0, size_type n1, size_type n2, RESIZE resize_allowed)
{
    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>
template<int R_,typename>
RA_INLINE_
void ra::rarray<T,R>::reshape(size_type n0,size_type n1,size_type n2,size_type n3, RESIZE resize_allowed)
{
    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>
template<int R_,typename>
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)
{
    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>
template<int R_,typename>
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)
{
    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>
template<int R_,typename>
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)
{
    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>
template<int R_,typename>
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)
{
    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>
template<int R_,typename>
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)
{
    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>
template<int R_,typename>
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)
{
    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>
template<int R_,typename>
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)
{
    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_NOEXCEPT(std::is_nothrow_copy_constructible<T>())
{
    RA_CHECKORSAY(not empty(), "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) RA_NOEXCEPT(true)
: 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>::empty() const noexcept
{
    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());
    return clone;
}


template<typename T, int R>
template<class U>
RA_INLINEF
typename std::enable_if<R!=1,U>::type
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>
template<class U>
RA_INLINEF
typename std::enable_if<R!=1,U>::type
ra::rarray<T,R>::at(size_type i) const
{
    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>
template<class U>
RA_INLINEF
typename std::enable_if<R==1,U>::type
ra::rarray<T,R>::at(size_type i)
{
    if (i < 0 or i >= extent(0))
        throw std::out_of_range("rarray<T,R>::at");
    return shape_.ptrs()[i];
}

template<typename T, int R>
template<class U>
RA_INLINEF
typename std::enable_if<R==1,U>::type
ra::rarray<T,R>::at(size_type i) const
{
    if (i < 0 or i >= extent(0))
        throw std::out_of_range("rarray<T,R>::at");
    return shape_.ptrs()[i];
}

template<typename T, int R> RA_INLINEF
typename ra::PointerArray<T,R-1>::type ra::rarray<T,R>::_at(ssize_t index)
{
    return shape_.ptrs()[index];
}

template<typename T, int R> RA_INLINEF
typename ra::PointerArray<T,R-1>::type ra::rarray<T,R>::_at(ssize_t index) const
{
    return shape_.ptrs()[index];
}


template<typename T, int R> RA_INLINEF
const ra::size_type* ra::rarray<T,R>::shape() const noexcept
{
    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 noexcept
{
    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 noexcept
{
    return const_cast<noconst_parray_t>(shape_.ptrs());
}

template<typename T, int R> RA_INLINE_
const typename ra::rarray<const T,R>& ra::rarray<T,R>::const_ref() const noexcept
{
    return reinterpret_cast<const 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) noexcept
{
    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) noexcept
{
    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 noexcept {
    return shape_.size();
}

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

template<typename T, int R> RA_INLINEF
const T* ra::rarray<T,R>::data() const noexcept
{
    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() noexcept
{
    return buffer_.begin();
}
template<typename T, int R> RA_INLINEF
typename ra::rarray<T,R>::const_iterator ra::rarray<T,R>::begin() const noexcept
{
    return buffer_.begin();
}
template<typename T, int R> RA_INLINEF
typename ra::rarray<T,R>::const_iterator ra::rarray<T,R>::cbegin() const noexcept
{
    return buffer_.cbegin();
}
template<typename T, int R> RA_INLINEF
typename ra::rarray<T,R>::iterator ra::rarray<T,R>::end() noexcept
{
    return buffer_.end();
}
template<typename T, int R> RA_INLINEF
typename ra::rarray<T,R>::const_iterator ra::rarray<T,R>::end() const noexcept
{
    return buffer_.end();
}
template<typename T, int R> RA_INLINEF
typename ra::rarray<T,R>::const_iterator ra::rarray<T,R>::cend() const noexcept
{
    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 T, int R, typename P> RA_INLINEF
ra::_Bracket<T,R-1,ra::_Bracket<T,R,P>> ra::_Bracket<T,R,P>::operator[](ssize_t nextindex) noexcept(noboundscheck)
{
    return { *this, nextindex, shape_ + 1 };
}
template<typename T, int R, typename P> RA_INLINEF
ra::_Bracket<T,R,P>::operator decltype(parent_.at(index_)) ()
{
    return parent_.at(index_);
}
template<typename T, int R, typename P> RA_INLINEF
ra::_Bracket<T,R,P>::_Bracket(P& parent, ssize_t index, const ssize_t* shape) noexcept(noboundscheck)
 : parent_(parent),
   index_(index),
   shape_(shape)
{
    RA_CHECKORSAY(index >=0 and index_ < shape_[0], "index out of range of array");
}
template<typename T, int R, typename P> RA_INLINEF
auto ra::_Bracket<T,R,P>::_at(ssize_t nextindex) noexcept(noboundscheck) -> decltype(parent_._at(index_)[nextindex])
{
    return parent_._at(index_)[nextindex];
}
template<typename T, int R, typename P> RA_INLINEF
auto ra::_Bracket<T,R,P>::at(ssize_t nextindex) -> decltype(parent_.at(index_).at(nextindex))
{
    return parent_.at(index_).at(nextindex);
}


template<typename T, typename P> RA_INLINEF
T& ra::_Bracket<T,1,P>::operator[](ssize_t nextindex) noexcept(noboundscheck)
{
    RA_CHECKORSAY(nextindex >=0 and nextindex < shape_[1], "index out of range of array");
    return parent_._at(index_)[nextindex];
}
template<typename T, typename P> RA_INLINEF
ra::_Bracket<T,1,P>::operator decltype(parent_.at(index_)) ()
{
    return parent_.at(index_);
}
template<typename T, typename P> RA_INLINEF
ra::_Bracket<T,1,P>::_Bracket(P& parent, ssize_t index, const ssize_t* shape) noexcept(noboundscheck)
 : parent_(parent),
   index_(index),
   shape_(shape)
{
    RA_CHECKORSAY(index >=0 and index_ < shape_[0], "index out of range of array");
}
template<typename T, typename P> RA_INLINEF
auto ra::_Bracket<T,1,P>::at(ssize_t nextindex) -> decltype(parent_.at(index_).at(nextindex))
{
    return parent_.at(index_).at(nextindex);
}


template<typename T,int R,typename P> RA_INLINEF
ra::_ConstBracket<T,R-1,ra::_ConstBracket<T,R,P>> ra::_ConstBracket<T,R,P>::operator[](ssize_t nextindex) const noexcept(noboundscheck)
{
    return { *this, nextindex, shape_ + 1 };
}
template<typename T, int R, typename P> RA_INLINEF
ra::_ConstBracket<T,R,P>::operator decltype(parent_.at(index_)) ()
{
    return parent_.at(index_);
}
template<typename T,int R,typename P> RA_INLINEF
ra::_ConstBracket<T,R,P>::_ConstBracket(const P& parent, ssize_t index, const ssize_t* shape) noexcept(noboundscheck)
: parent_(parent),
  index_(index),
  shape_(shape)
{
    RA_CHECKORSAY(index >=0 and index_ < shape_[0], "index out of range of array");
}
template<typename T,int R,typename P> RA_INLINEF
auto ra::_ConstBracket<T,R,P>::_at(ssize_t nextindex) const noexcept(noboundscheck) -> decltype(parent_._at(index_)[nextindex])
{
    return parent_._at(index_)[nextindex];
}
template<typename T, int R, typename P> RA_INLINEF
auto ra::_ConstBracket<T,R,P>::at(ssize_t nextindex) const -> decltype(parent_.at(index_).at(nextindex))
{
    return parent_.at(index_).at(nextindex);
}


template<typename T,typename P> RA_INLINEF
const T& ra::_ConstBracket<T,1,P>::operator[](ssize_t nextindex) const noexcept(noboundscheck)
{
    RA_CHECKORSAY(nextindex >=0 and nextindex < shape_[1], "index out of range of array");
    return parent_._at(index_)[nextindex];
}
template<typename T, typename P> RA_INLINEF
ra::_ConstBracket<T,1,P>::operator decltype(parent_.at(index_)) ()
{
    return parent_.at(index_);
}
template<typename T,typename P> RA_INLINEF
ra::_ConstBracket<T,1,P>::_ConstBracket(const P& parent, ssize_t index, const ssize_t* shape) noexcept(noboundscheck)
: parent_(parent),
  index_(index),
  shape_(shape)
{
    RA_CHECKORSAY(index >=0 and index_ < shape_[0], "index out of range of array");
}
template<typename T, typename P> RA_INLINEF
auto ra::_ConstBracket<T,1,P>::at(ssize_t nextindex) const -> decltype(parent_.at(index_).at(nextindex))
{
    return parent_.at(index_).at(nextindex);
}


namespace ra {
template<typename A>
size_type extent(const A &a, int i)
{
    switch (i) {
    case 0: return std::extent<A,0>();
    case 1: return std::extent<A,1>();
    case 2: return std::extent<A,2>();
    case 3: return std::extent<A,3>();
    case 4: return std::extent<A,4>();
    case 5: return std::extent<A,5>();
    case 6: return std::extent<A,6>();
    case 7: return std::extent<A,7>();
    case 8: return std::extent<A,8>();
    case 9: return std::extent<A,9>();
    case 10: return std::extent<A,10>();
    default: return 0;
    }
}

template<typename T, int R>
size_type extent(const rarray<T,R> &a, int i)
{
    return a.extent(i);
}
}


//begin #include "rarrayio.h"

namespace ra {

template<typename T,int R> RA_INLINE_ std::ostream& text_output(std::ostream &o, const rarray<T,R>& r);
template<typename T>       RA_INLINE_ std::ostream& text_output(std::ostream &o, const rarray<T,1>& r);

template<typename T, int R>
struct Deref {
    static RA_INLINEF T& access(typename PointerArray<T,R>::type p, const size_type* indices);
};
template<typename T>
struct Deref<T,1>
{
    static RA_INLINEF T& access(typename PointerArray<T,1>::type p, const size_type* indices);
};

template<typename T>
struct StringToValue {
    static RA_INLINE_ void get(const std::string& input, T& output);
};
template<>
struct StringToValue<std::string> {
    static RA_INLINE_ void get(const std::string& input, std::string& output);
};

enum class token { UNKNOWN=0, BRACEOPEN, BRACECLOSE, COMMA, DATASTRING, END };

 inline char toch(const token& Token) {
     switch (Token) {
     case ra::token::BRACEOPEN:  return '{';
     case ra::token::BRACECLOSE: return '}';
     case ra::token::COMMA:      return ',';
     case ra::token::DATASTRING: return '$';
     case ra::token::END:        return '.';
     default:                    return '?';
     }
 }

RA_INLINE_
std::list<std::pair<token,std::string>> parse_shape(std::istream & in, int R, size_type* shape);

template<typename T, int R> RA_INLINE_
void parse_strings(const std::list<std::pair<token,std::string>> & tokens, const size_type* shape, typename PointerArray<T,R>::type p);

}



template<typename T, int R> RA_INLINE_ std::ostream& ra::operator<<(std::ostream &o, const ra::rarray<T, R>& r)
{
    if (R>1) {
        return ra::text_output(o,r);
    } else if (R==1) {
        if (not r.empty()) {
            o << "{";
            for (ra::size_type i=0; i<r.extent(0); i++)  {
                if (i>0)
                    o << ",";
                o << r.data()[i];
            }
            o << "}";
        }
    }
    return o;
}


template<typename T,int R> RA_INLINE_
std::ostream& ra::text_output(std::ostream &o, const ra::rarray<T,R>& r)
{
    if (not r.empty()) {
        o << "{\n";
        for (size_type i=0; i<r.extent(0); i++)  {
            if (i>0)
                o << ",\n";
            o << r.at(i);
        }
        o << "\n}";
    } else {
        for (int i=0; i<R; i++)
            o << '{';
        for (int i=0; i<R; i++)
            o << '}';
        o << "\n";
    }
    return o;
}

template<typename T> RA_INLINE_
std::ostream& ra::text_output(std::ostream &o, const ra::rarray<T,1>& r)
{
    if (not r.empty()) {
        o << '{';
        for (size_type i=0; i<r.extent(0); i++) {
            if (i) o << ',';
            std::stringstream strstr;
            std::string result;
            const T& val = r.at(i);
            strstr << val;
            result = strstr.str();
            if (result.find_first_of("{,}#") != std::string::npos
                and not
                (result[0]=='(' and result[result.size()-1]==')' and result.substr(1,result.size()-2).find_first_of(")") == std::string::npos) )
                o << '#' << result.size() << ':';
            o << result;
        }
        o << '}';
    } else
        o << "{}";
    return o;
}



template<typename T> RA_INLINE_
void ra::StringToValue<T>::get(const std::string& input, T& output)
{
    std::stringstream str(input);
    str >> output;
}

RA_INLINE_ void ra::StringToValue<std::string>::get(const std::string& input, std::string& output)
{
    output = input;
}


template<typename T, int R> RA_INLINEF
T& ra::Deref<T,R>::access(typename PointerArray<T,R>::type p, const size_type* indices)
{
    return Deref<T,R-1>::access(p[indices[0]-1], indices+1);
}

template<typename T> RA_INLINEF
T& ra::Deref<T,1>::access(typename PointerArray<T,1>::type p, const size_type* indices)
{
    return p[indices[0]-1];
}


static char get_but_eat_newline(std::istream & in)
{
    char ch1='\n';
    while (ch1=='\n' && ! in.eof())
        in >> ch1;
    return ch1;
}

static char get_but_eat_whitespace(std::istream & in)
{
    char ch1;
    in >> ch1;
    return ch1;
}

RA_INLINE_
std::list<std::pair<ra::token,std::string>> ra::parse_shape(std::istream & in, int R, ra::size_type* shape)
{
    if (shape == nullptr) {
        throw std::istream::failure("Parse call error");
    }
    std::list<std::pair<ra::token,std::string>> result;
    size_t init_file_ptr = in.tellg();
    try {
        size_type current_shape[R];
        for (int i=0; i<R; i++) {
            current_shape[i] = 1;
            shape[i] = 0;
            if (get_but_eat_newline(in) != '{')
                throw std::istream::failure("Format error");
            else
                result.push_back({ra::token::BRACEOPEN,""});
        }
        int current_depth = R-1;
        while (current_depth>=0) {
            if (current_depth==R-1) {
                char         lastchar;
                std::string  word = "";
                do {
                    if (word=="") {
                        lastchar = get_but_eat_newline(in);
                    } else {
                        in.get(lastchar);
                    }
                    if (lastchar != ',' and lastchar != '}') {
                        word += lastchar;
                    }
                    if (word == "#") {
                        word="";
                        std::string skipstr;
                        do {
                            in.get(lastchar);
                            skipstr += lastchar;
                        } while (lastchar!=':');
                        int skip = atoi(skipstr.c_str());//
                        for (int i=0; i<skip; i++) { //
                            char nextchar;
                            in.get(nextchar);
                            word += nextchar;
                        }
                    } else if (word == "(") {
                        const int safeguardcount=1024*1024; //
                        int count=0; //
                        while (lastchar!=')' and count<safeguardcount) {
                            in.get(lastchar);
                            word += lastchar;
                            count++;
                        }
                        in.get(lastchar);
                    }
                    if (lastchar == ',') {
                        result.push_back({ra::token::DATASTRING,word});
                        result.push_back({ra::token::COMMA,""});
                        word="";
                        current_shape[current_depth]++;
                    }
                } while (lastchar != '}');
                result.push_back({ra::token::DATASTRING,word});
                result.push_back({ra::token::BRACECLOSE,""});
                if (shape)
                    if (shape[current_depth] < current_shape[current_depth])
                        shape[current_depth] = current_shape[current_depth];
                current_depth--;
            } else {
                switch (get_but_eat_whitespace(in)) {
                   case ',':
                    result.push_back({ra::token::COMMA,""});
                    current_shape[current_depth]++;
                    break;
                  case '{':
                    result.push_back({ra::token::BRACEOPEN,""});
                    current_depth++;
                    current_shape[current_depth] = 1;
                    break;
                  case '}':
                    result.push_back({ra::token::BRACECLOSE,""});
                    if (shape)
                          if (shape[current_depth] < current_shape[current_depth])
                              shape[current_depth] = current_shape[current_depth];
                    current_depth--;
                    break;
                  default:
                    throw std::istream::failure("Format error");
                }
            }
        }
        result.push_back({ra::token::END,""});
    }
    catch (std::istream::failure& e) {
        in.seekg(init_file_ptr, in.beg);
        in.setstate(std::ios::failbit);
    }

    return result;
}

template<typename T, int R> RA_INLINE_
void ra::parse_strings(const std::list<std::pair<token,std::string>> & tokens, const ra::size_type* shape, typename ra::PointerArray<T,R>::type p)
{
    size_type index[R];
    int current_depth = -1;
    for (auto& tokenpair: tokens) {
        switch (tokenpair.first) {
        case ra::token::BRACEOPEN:
            current_depth++;
            index[current_depth]=1;
            break;
        case ra::token::BRACECLOSE:
            current_depth--;
            break;
        case ra::token::COMMA:
            index[current_depth]++;
            break;
        case ra::token::DATASTRING:
            StringToValue<T>::get(tokenpair.second, Deref<T,R>::access(p, index));
            break;
        case ra::token::END:
            break;
        default:
            throw std::istream::failure("Parsing error");
            break;
        }
        if (tokenpair.first == ra::token::END) break;
    }
}

template<typename T,int R> RA_INLINE_
std::istream& ra::operator>>(std::istream &in, ra::rarray<T,R>& r)
{
    ra::size_type extent[R] = {0};
    auto X = ra::parse_shape(in, R, extent);
    if (ra::mul(extent,R) <= r.size())
        r.reshape(extent, ra::RESIZE::ALLOWED);
    else
        r = rarray<T,R>(extent);
    ra::parse_strings<T,R>(X, extent, r.ptr_array());
    return in;
}


//end #include "rarrayio.h"


//begin #include "rarraydelmacros.h"

#undef RA_CHECKORSAY
#undef RA_NOEXCEPT
#undef RA_INLINE
#undef RA_INLINEF
#undef RA_INLINE_
#undef noboundscheck


//end #include "rarraydelmacros.h"



#define EXTENT(A,I)  ra::extent(A,I)
#define RARRAY(A)    rarray<std::remove_all_extents<decltype(A)>::type,std::rank<decltype(A)>::value>(A)
#define INDEX(A,X,I) RARRAY(A).index(X,I)

namespace std {
    template<class T, int R >
    struct remove_all_extents<ra::rarray<T,R>> {
        typedef T type;
    };
    template<class T, int R >
    struct rank<ra::rarray<T,R>> {
        static const size_t value = R;
    };
}

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