root/branches/newmole/source/container_classes.h

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */

#ifndef _CONTAINER_CLASSES_H_
#define _CONTAINER_CLASSES_H_

void do_dump_state(const void* buf, size_t nelem, size_t size, FILE* out, int32 magic);
void do_restore_state(void* buf, size_t nelem, size_t size, FILE *in, int32 magic);

//! define the supported memory layouts of multi_arr / flex_arr
//!
//! ARPA_TYPE: arrays of pointers to arrays to ...
//! C_TYPE: row-major order, exactly as used in C (last index runs fastest)
//! FLX_TYPE: the one and only layout used by flex_arr
typedef enum { ARPA_TYPE, C_TYPE, FLX_TYPE, ML_TOP } mem_layout;

// magic numbers to identify each memory layout
static const int32 MA_VERS[ML_TOP] = { 120070905, 220070803, 320071126 };

#ifdef USE_C_TYPE
#define MEM_LAYOUT_VAL C_TYPE
#else
#define MEM_LAYOUT_VAL ARPA_TYPE
#endif

#ifdef BOUNDS_CHECK
#define lgBOUNDSCHECKVAL true
#else
#define lgBOUNDSCHECKVAL false
#endif

//! basic_pntr - base class for generalization of normal pointers that enables bounds checking
//! it comes with the full set of operators that you would expect for a random access pointer
template<class T, bool lgBC>
class basic_pntr
{
        static const int p_nd = lgBC ? 3 : 1;
        T* p_p[p_nd]; // p[0] current pointer, p[1] lower bound, p[2] upper bound

        T* p_index_checked ( const ptrdiff_t n ) const
        {
                T* t = p_p[0]+n;
#ifdef _MSC_VER
                /* disable warning that conditional expression is constant, true or false in if */
#               pragma warning( disable : 4127 )
#endif
                if( lgBC ) 
                {
                        if( t < p_p[1] || t >= p_p[2] )
                                OUT_OF_RANGE( "basic_pntr::p_index_checked()" );
                }
                return t;
        }
        void p_set_vals( T* p0, T* p1, T* p2 )
        {
#ifdef _MSC_VER
                /* disable warning that conditional expression is constant, true or false in if */
#               pragma warning( disable : 4127 )
#endif
                if( lgBC )
                {
                        p_p[0] = p0; p_p[1] = p1; p_p[2] = p2;
                }
                else
                        p_p[0] = p0;
        }
public:
        typedef random_access_iterator_tag iterator_category;
        typedef T            value_type;
        typedef T&           reference;
        typedef const T&     const_reference;
        typedef T*           pointer;
        typedef const T*     const_pointer;
        typedef size_t       size_type;
        typedef ptrdiff_t    difference_type;

        // constructors
        basic_pntr( T* p0, T* p1, T* p2 )
        {
                p_set_vals( p0, p1, p2 );
        }
        basic_pntr( T* p0 )
        {
                p_set_vals( p0, NULL, NULL );
        }
        basic_pntr()
        {
                p_set_vals( NULL, NULL, NULL );
        }
        basic_pntr( const basic_pntr& t )
        {
                *this = t;
        }
        // protected destructor (to prevent direct instantiation or destruction via basic_pntr*, see Alexandrescu p13)
protected:
        ~basic_pntr() {} 
public:
        // pre-increment
        basic_pntr& operator++ ()
        {
                ++p_p[0];
                return *this;
        }
        // pre-decrement
        basic_pntr& operator-- ()
        {
                --p_p[0];
                return *this;
        }
        // define operators for += and -=, normal arithmetic is defined separately in
        // the derived classes; it cannot be done here since they would require implicit
        // conversion from basic_pntr -> pntr or const_pntr to work; this would also create
        // an implicit and silent conversion from const_pntr -> pntr, which is illegal...
        basic_pntr& operator+= ( const ptrdiff_t n ) { p_p[0] += n; return *this; }
        basic_pntr& operator-= ( const ptrdiff_t n ) { p_p[0] -= n; return *this; }
        // dereference
        T& operator* () const
        {
                return *(p_index_checked(0));
        }
        T* operator-> () const
        {
                return p_index_checked(0);
        }
        T& operator[] ( const ptrdiff_t n ) const
        {
                return *(p_index_checked(n));
        }
        // finally, define the boolean operators...
        bool operator== ( const basic_pntr& t ) const { return p_p[0] == t.p_p[0]; }
        bool operator!= ( const basic_pntr& t ) const { return p_p[0] != t.p_p[0]; }
        bool operator<  ( const basic_pntr& t ) const { return p_p[0] <  t.p_p[0]; }
        bool operator<= ( const basic_pntr& t ) const { return p_p[0] <= t.p_p[0]; }
        bool operator>  ( const basic_pntr& t ) const { return p_p[0] >  t.p_p[0]; }
        bool operator>= ( const basic_pntr& t ) const { return p_p[0] >= t.p_p[0]; }
};

#ifdef __SUNPRO_CC
// these beauties are needed by the Rogue Wave implementation of STL algorithms used by the
// Sun Studio compiler, otherwise their STL algorithms won't compile with pntr or const_pntr
namespace std
{
        template<class T, bool lgBC>
        inline T* __value_type(const basic_pntr<T,lgBC>&)
        {
                return static_cast<T*>( NULL );
        }

        template<class T, bool lgBC>
        inline typename basic_pntr<T,lgBC>::difference_type* __distance_type(const basic_pntr<T,lgBC>&)
        {
                return static_cast<basic_pntr<T,lgBC>::difference_type*>( NULL );
        }
}
#endif

//! pntr - interface class to replace normal pointers
template<class T, bool lgBC>
class pntr : public basic_pntr<T,lgBC>
{
public:
        // constructors are not inherited, so define them again
        pntr( T* p0 ) : basic_pntr<T,lgBC>( p0 ) {}
        pntr( T* p0, T* p1, T* p2 ) : basic_pntr<T,lgBC>( p0, p1, p2 ) {}
        pntr() {}
        // the increment / decrement operators need to be recast...
        // otherwise expressions like p = ++q would be illegal for iterators...
        pntr& operator++ () { return static_cast<pntr&>(basic_pntr<T,lgBC>::operator++()); }
        const pntr operator++ (int) { pntr t = *this; ++(*this); return t; }
        pntr& operator-- () { return static_cast<pntr&>(basic_pntr<T,lgBC>::operator--()); }
        const pntr operator-- (int) { pntr t = *this; --(*this); return t; }
        // define p+n, p-n, p-q
        const pntr operator+ ( const ptrdiff_t n ) const { pntr s = *this; s += n; return s; }
        const pntr operator- ( const ptrdiff_t n ) const { pntr s = *this; s -= n; return s; }
        ptrdiff_t operator- ( const pntr& t ) const { return &(*this[0]) - &t[0]; }
};

// this defines n+p
template<class T, bool lgBC>
inline const pntr<T,lgBC> operator+ ( const ptrdiff_t n, const pntr<T,lgBC>& t )
{
        pntr<T,lgBC> s = t;
        s += n;
        return s;
}

//! const_pntr - same as pntr, except that it replaces const pointers rather than normal pointers
template<class T, bool lgBC>
class const_pntr : public basic_pntr<T,lgBC>
{
public:
        // constructors are not inherited, so define them again
        const_pntr( T* p0 ) : basic_pntr<T,lgBC>( p0 ) {}
        const_pntr( T* p0, T* p1, T* p2 ) : basic_pntr<T,lgBC>( p0, p1, p2 ) {}
        const_pntr() {}
        // make sure we can assign a pntr to a const_pntr by creating an implicit conversion to const_pntr
        const_pntr( const pntr<T,lgBC>& t ) : basic_pntr<T,lgBC>( t ) {}
        // the increment / decrement operators need to be recast...
        // otherwise expressions like *p++ = 1. would be legal for const_iterators...
        const_pntr& operator++ () { return static_cast<const_pntr&>(basic_pntr<T,lgBC>::operator++()); }
        const const_pntr operator++ (int) { const_pntr t = *this; ++(*this); return t; }
        const_pntr& operator-- () { return static_cast<const_pntr&>(basic_pntr<T,lgBC>::operator--()); }
        const const_pntr operator-- (int) { const_pntr t = *this; --(*this); return t; }
        const_pntr& operator+= ( const ptrdiff_t n )
        {
                return static_cast<const_pntr&>(basic_pntr<T,lgBC>::operator+=(n));
        }
        const_pntr& operator-= ( const ptrdiff_t n )
        {
                return static_cast<const_pntr&>(basic_pntr<T,lgBC>::operator-=(n));
        }
        // the dereference operators need to be recast...
        const T& operator* () const { return static_cast<const T&>(basic_pntr<T,lgBC>::operator*()); }
        const T* operator-> () const { return static_cast<const T*>(basic_pntr<T,lgBC>::operator->()); }
        const T& operator[] ( const ptrdiff_t n ) const
        {
                return static_cast<const T&>(basic_pntr<T,lgBC>::operator[](n));
        }
        // define p+n, p-n, p-q
        const const_pntr operator+ ( const ptrdiff_t n ) const { const_pntr s = *this; s += n; return s; }
        const const_pntr operator- ( const ptrdiff_t n ) const { const_pntr s = *this; s -= n; return s; }
        ptrdiff_t operator- ( const const_pntr& t ) const { return &(*this[0]) - &t[0]; }
};

// this defines n+p
template<class T, bool lgBC>
inline const const_pntr<T,lgBC> operator+ ( const ptrdiff_t n, const const_pntr<T,lgBC>& t )
{
        const_pntr<T,lgBC> s = t;
        s += n;
        return s;
}

//! tree_vec - a simple class to store the bounds checking information for multi_arr
struct tree_vec
{
        typedef size_t size_type;

        size_type n;
        tree_vec *d;

private:
        void p_clear0()
        {
                if( d != NULL )
                {
                        for( size_type i = 0; i < n; ++i )
                                d[i].clear();
                        delete[] d;
                }
        }
        void p_clear1()
        {
                n = 0;
                d = NULL;
        }

public:
        tree_vec()
        {
                p_clear1();
        }
        ~tree_vec()
        {
                p_clear0();
        }
        void clear()
        {
                p_clear0();
                p_clear1();
        }
        const tree_vec& operator= ( const tree_vec& m )
        {
                if( &m != this )
                {
                        clear();
                        n = m.n;
                        if( m.d != NULL )
                        {
                                d = new tree_vec[n];
                                tree_vec *p = d;
                                const tree_vec *mp = m.d;
                                for( size_type i = 0; i < n; ++i )
                                        *p++ = *mp++;
                        }
                }
                return *this;
        }
        tree_vec& getvec(const size_type i, const size_type index[])
        {
                if( i == 0 )
                        return *this;
                else
                        return getvec(i-1,index).d[index[i-1]];
        }
        const tree_vec& getvec(const size_type i, const size_type index[]) const
        {
                if( i == 0 )
                        return *this;
                else
                        return getvec(i-1,index).d[index[i-1]];
        }
};

//! multi_geom - this class maintains all the geometry information for multi_arr
//! keeping it separate makes it easy to clone the information from one multi_arr to another
template<int d,mem_layout ALLOC=MEM_LAYOUT_VAL>
class multi_geom
{
public:
        typedef size_t size_type;

        tree_vec v;

        size_type size;  //! allocated size (number of data elements, pointers are not counted)
        size_type s[d];  //! size of each dimension (only used in C_TYPE layout)
        size_type st[d]; //! stride for each dimension (only used in C_TYPE layout)
        size_type nsl[d];//! sizes of each of the pointer arrays

private:
        void p_clear0()
        {
                v.clear();
        }
        void p_clear1()
        {
                size = 0;
                for( int i=0; i < d; ++i )
                {
                        s[i] = 0;
                        st[i] = 0;
                        nsl[i] = 0;
                }
        }

public:
        multi_geom()
        {
                p_clear1();
        }
        ~multi_geom()
        {
                p_clear0();
        }
        void clear()
        {
                p_clear0();
                p_clear1();
        }
        const multi_geom& operator= ( const multi_geom& m )
        {
                if( &m != this )
                {
                        clear();
                        v = m.v;
                        size = m.size;
                        for( int i=0; i < d; ++i )
                        {
                                s[i] = m.s[i];
                                st[i] = m.st[i];
                                nsl[i] = m.nsl[i];
                        }
                }
                return *this;
        }
        bool lgInbounds(const size_type n, const size_type index[]) const
        {
                if( n != 0 )
                        return ( lgInbounds(n-1,index) && index[n-1] < v.getvec(n-1,index).n );
                else
                        return true;
        }
        void reserve(const size_type n, const size_type index[])
        {
                ASSERT( n <= d && index[n-1] > 0 && lgInbounds( n-1, index ) );

                tree_vec& w = v.getvec( n-1, index );
                if( d > n )
                {
                        ASSERT( w.d == NULL );
                        w.d = new tree_vec[ index[n-1] ];
                }
                w.n = index[n-1];
                s[n-1] = max(s[n-1],index[n-1]);
                nsl[n-1] += index[n-1];
        }
        void reserve_recursive(const size_type n, size_type index[])
        {
                if( n == 0 )
                {
                        reserve( n+1, index );
                        if( n+1 < d )
                                reserve_recursive( n+1, index );
                }
                else
                {
                        size_type top = index[n-1];
                        for( size_type i=0; i < top; ++i )
                        {
                                index[n-1] = i;
                                reserve( n+1, index );
                                if( n+1 < d )
                                        reserve_recursive( n+1, index );
                        }
                        index[n-1] = top;
                }
        }
        void finalize(void)
        {
#ifdef _MSC_VER
                /* disable warning that conditional expression is constant, true or false in if */
#       pragma warning( disable : 4127 )
#endif

                STATIC_ASSERT ( ALLOC == C_TYPE || ALLOC == ARPA_TYPE );

                if( ALLOC == ARPA_TYPE )
                {
                        size_type n1[d], n2[d];
                        for( int dim=0; dim < d; ++dim )
                                n1[dim] = n2[dim] = 0L;
                        // sanity checks
                        p_setupArray( n1, n2, &v, 0 );
                        for( int dim=0; dim < d-1; ++dim )
                                ASSERT( n1[dim] == nsl[dim] && n2[dim] == nsl[dim+1] );
                        size = nsl[d-1];
                }
                else if( ALLOC == C_TYPE )
                {
                        st[d-1] = s[d-1]; 
                        for( int i = d-2; i >= 0; --i )
                                st[i] = st[i+1]*s[i];
                        size = st[0];
                }
                else
                {
                        TotalInsanity();
                }
        }

private:
        void p_setupArray( size_type n1[], size_type n2[], const tree_vec* w, int l )
        {
                for( size_type i=0; i < w->n; ++i )
                {
                        n1[l]++;
                        if( l < d-2 )
                        {
                                p_setupArray( n1, n2, &w->d[i], l+1 );
                        }
                        n2[l] += w->d[i].n;
                }
        }
};


//
//! The n_pointer and const_n_pointer classes below define indexing into the multi_arr
//!
//! NB NB -- The design of these classes is CRUCIAL for CPU performance! Any redesign
//!          should be thoroughly tested for CPU time regressions on a broad range of
//!          platforms and compilers!
//!
//! The primary goal of the design of these classes should be to make life as easy as
//! possible for the compiler, NOT the programmer! The design should assure that that
//! the methods are FULLY INLINED. Without that the resulting code will be severely
//! crippled!
//!
//! Below are a total of 8 specializations for each of these classes. These could be
//! combined into 2 specialization for each class (one for N > 1 and one for N == 1),
//! but it turns out that the resulting method is too complicated for most compilers
//! to inline. Hence we decided to split the methods up into simpler specializations.
//! As they are written below, nearly every compiler can inline them (at least for the
//! non-bounds_checking versions)...
//!
//! This is the situation on 2008 jan 26:
//!
//! Linux IA32/AMD64   ARPA    C    ARPA/BC  C/BC
//!    g++ 4.3.x        OK     OK     OK      OK   (tested on prerelease)
//!    g++ 4.2.x        OK     OK     OK      OK
//!    g++ 4.1.x        OK     OK     OK      OK
//!    g++ 4.0.x        OK     OK     OK      OK
//!    g++ 3.4.x        ??     OK     ??      ??
//!    g++ 3.3.x        ??     OK     ??      ??
//!
//!    icc 10.1         OK     OK     OK      OK
//!    icc 10.0         OK     OK     OK      OK
//!    icc  9.1         OK     OK     OK      OK
//!    icc  9.0         OK     OK     OK      OK
//!
//!    Sun Studio 12    OK     OK     OK      OK
//!
//!    Pathscale 2.3.1  OK     OK     --      --   (only tested on AMD64)
//!
//!    Portland 6.1-5   --     --     --      --   (only tested on IA32)
//!
//! Linux IA64
//!    g++ 4.2.x        OK     OK     OK      OK
//!    g++ 3.3.x        ??     OK?    ??      ??
//!
//!    icc 10.1         OK     OK     --      --
//!    icc 10.0         OK     OK     --      --
//!    icc  9.1         OK     OK     --      --
//!    icc  9.0         OK     OK     OK      OK
//!    icc  8.1         OK     OK     OK      OK
//!
//! Windows IA32
//!    Visual Studio 8  OK     OK     --      --
//!    icl 10.1         OK     OK     --      --
//!
//! Solaris Ultrasparc
//!    g++ 4.3.x        OK     OK     OK      OK   (tested on prerelease)
//!    g++ 4.2.x        OK     OK     OK      OK
//!    g++ 4.1.x        OK     OK     OK      OK
//!    g++ 4.0.x        OK     OK     OK      OK
//!    g++ 3.4.x        ??     ??     ??      ??
//!    g++ 3.3.x        ??     ??     ??      ??
//!
//!    Sun Studio 12    OK     OK     OK      OK
//!
//! Notes:
//!   - The tests were carried out on the following simple test code:
//!
//!       long test(multi_arr<long,6>& arr)
//!       {
//!         return arr[2][3][4][5][6][7];
//!       }
//!
//!     The tests were done both on ARPA_TYPE and C_TYPE arrays, with and without
//!     bounds-checking enabled (BC indicates bounds-checking enabled). The resulting
//!     (optimized) assembly was inspected by eye.
//!
//!   - "OK" indicates that all methods were inlined and the resulting assembly looked
//!     optimal; "??" indicates that all methods were inlined, but the assembly looked
//!     sub-optimal (i.e. not all overhead was optimized away); "--" indicates that at
//!     least some of the methods were not inlined at all.
//!
//! Indexing a multi_arr works as follows. When the compiler encounters arr[i][j][k],
//! it will first emit a to call multi_arr::operator[] with i as its argument; this
//! operator will construct a temporary n_pointer as follows
//!
//!     arr[i]  -->  n_pointer<T,3>(*p,st[],*v)::operator[](i) -> n_pointer<T,2>
//!
//! here p is the base pointer to the data, st contains the strides of a C_TYPE array
//! (which are not used in ARPA_TYPE arrays) and v points to the tree_vec with the
//! bounds-checking information. The operator[] of the n_pointer will update p and v
//! using the value of i. This will only take care of the first index. At this stage
//! the compiler still needs to emit code for the remainder: n_pointer<T,2>[j][k].
//! The operator[] of the n_pointer<T,2> will take care of the next index:
//!
//!     n_pointer<T,2>::operator[](j)  -->  n_pointer<T,1>
//!
//! So another temporary n_pointer is emitted, and p and v are updated again by the
//! operator[] using j. Now the compiler still needs to emit code for: n_pointer<T,1>[k];
//! n_pointer<T,1> is special since it will absorb the last index, so it should not
//! return yet another n_pointer but a (const) reference to the data item itself:
//!
//!     n_pointer<T,1>::operator[](k)  -->  T& *(p + k)
//!     const_n_pointer<T,1>::operator[](k)  -->  const T& *(p + k)
//


// forward definitions
template<class T, int N, mem_layout ALLOC, bool lgBC> class n_pointer;
template<class T, int N, mem_layout ALLOC, bool lgBC> class const_n_pointer;

template<class T, int N>
class n_pointer<T,N,ARPA_TYPE,false>
{
        T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        n_pointer(T* p, const size_t* st=NULL, const tree_vec* v=NULL) : p_p(p), p_st(st), p_v(v) {}
        const n_pointer<T,N-1,ARPA_TYPE,false> operator[] (const size_t i) const
        {
                return n_pointer<T,N-1,ARPA_TYPE,false>( *((T**)p_p+i) );
        }
};

template<class T, int N>
class n_pointer<T,N,C_TYPE,false>
{
        T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        n_pointer(T* p, const size_t* st, const tree_vec* v=NULL) : p_p(p), p_st(st), p_v(v) {}
        const n_pointer<T,N-1,C_TYPE,false> operator[] (const size_t i) const
        {
                return n_pointer<T,N-1,C_TYPE,false>( p_p+i*p_st[0], p_st+1 );
        }
};

template<class T>
class n_pointer<T,1,ARPA_TYPE,false>
{
        T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        n_pointer(T* p, const size_t* st=NULL, const tree_vec* v=NULL) : p_p(p), p_st(st), p_v(v) {}
        T& operator[] (const size_t i) const
        {
                return *(p_p + i);
        }
};

template<class T>
class n_pointer<T,1,C_TYPE,false>
{
        T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        n_pointer(T* p, const size_t* st, const tree_vec* v=NULL) : p_p(p), p_st(st), p_v(v) {}
        T& operator[] (const size_t i) const
        {
                return *(p_p + i);
        }
};

template<class T, int N>
class n_pointer<T,N,ARPA_TYPE,true>
{
        T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        n_pointer(T* p, const size_t* st, const tree_vec* v) : p_p(p), p_st(st), p_v(v) {}
        const n_pointer<T,N-1,ARPA_TYPE,true> operator[] (const size_t i) const
        {
                if( i >= p_v->n )
                        OUT_OF_RANGE( "n_pointer::operator[]" );
                return n_pointer<T,N-1,ARPA_TYPE,true>( *((T**)p_p+i), NULL, &p_v->d[i] );
        }
};

template<class T, int N>
class n_pointer<T,N,C_TYPE,true>
{
        T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        n_pointer(T* p, const size_t* st, const tree_vec* v) : p_p(p), p_st(st), p_v(v) {}
        const n_pointer<T,N-1,C_TYPE,true> operator[] (const size_t i) const
        {
                if( i >= p_v->n )
                        OUT_OF_RANGE( "n_pointer::operator[]" );
                return n_pointer<T,N-1,C_TYPE,true>( p_p+i*p_st[0], p_st+1, &p_v->d[i] );
        }
};

template<class T>
class n_pointer<T,1,ARPA_TYPE,true>
{
        T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        n_pointer(T* p, const size_t* st, const tree_vec* v) : p_p(p), p_st(st), p_v(v) {}
        T& operator[] (const size_t i) const
        {
                if( i >= p_v->n )
                        OUT_OF_RANGE( "n_pointer::operator[]" );
                return *(p_p + i);
        }
};

template<class T>
class n_pointer<T,1,C_TYPE,true>
{
        T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        n_pointer(T* p, const size_t* st, const tree_vec* v) : p_p(p), p_st(st), p_v(v) {}
        T& operator[] (const size_t i) const
        {
                if( i >= p_v->n )
                        OUT_OF_RANGE( "n_pointer::operator[]" );
                return *(p_p + i);
        }
};

template<class T, int N>
class const_n_pointer<T,N,ARPA_TYPE,false>
{
        const T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        const_n_pointer(const T* p, const size_t* st=NULL, const tree_vec* v=NULL) : p_p(p), p_st(st), p_v(v) {}
        const const_n_pointer<T,N-1,ARPA_TYPE,false> operator[] (const size_t i) const
        {
                return const_n_pointer<T,N-1,ARPA_TYPE,false>( *((T**)p_p+i) );
        }
};

template<class T, int N>
class const_n_pointer<T,N,C_TYPE,false>
{
        const T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        const_n_pointer(const T* p, const size_t* st, const tree_vec* v=NULL) : p_p(p), p_st(st), p_v(v) {}
        const const_n_pointer<T,N-1,C_TYPE,false> operator[] (const size_t i) const
        {
                return const_n_pointer<T,N-1,C_TYPE,false>( p_p+i*p_st[0], p_st+1 );
        }
};

template<class T>
class const_n_pointer<T,1,ARPA_TYPE,false>
{
        const T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        const_n_pointer(const T* p, const size_t* st=NULL, const tree_vec* v=NULL) : p_p(p), p_st(st), p_v(v) {}
        const T& operator[] (const size_t i) const
        {
                return *(p_p + i);
        }
};

template<class T>
class const_n_pointer<T,1,C_TYPE,false>
{
        const T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        const_n_pointer(const T* p, const size_t* st, const tree_vec* v=NULL) : p_p(p), p_st(st), p_v(v) {}
        const T& operator[] (const size_t i) const
        {
                return *(p_p + i);
        }
};

template<class T, int N>
class const_n_pointer<T,N,ARPA_TYPE,true>
{
        const T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        const_n_pointer(const T* p, const size_t* st, const tree_vec* v) : p_p(p), p_st(st), p_v(v) {}
        const const_n_pointer<T,N-1,ARPA_TYPE,true> operator[] (const size_t i) const
        {
                if( i >= p_v->n )
                        OUT_OF_RANGE( "const_n_pointer::operator[]" );
                return const_n_pointer<T,N-1,ARPA_TYPE,true>( *((T**)p_p+i), NULL, &p_v->d[i] );
        }
};

template<class T, int N>
class const_n_pointer<T,N,C_TYPE,true>
{
        const T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        const_n_pointer(const T* p, const size_t* st, const tree_vec* v) : p_p(p), p_st(st), p_v(v) {}
        const const_n_pointer<T,N-1,C_TYPE,true> operator[] (const size_t i) const
        {
                if( i >= p_v->n )
                        OUT_OF_RANGE( "const_n_pointer::operator[]" );
                return const_n_pointer<T,N-1,C_TYPE,true>( p_p+i*p_st[0], p_st+1, &p_v->d[i] );
        }
};

template<class T>
class const_n_pointer<T,1,ARPA_TYPE,true>
{
        const T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        const_n_pointer(const T* p, const size_t* st, const tree_vec* v) : p_p(p), p_st(st), p_v(v) {}
        const T& operator[] (const size_t i) const
        {
                if( i >= p_v->n )
                        OUT_OF_RANGE( "const_n_pointer::operator[]" );
                return *(p_p + i);
        }
};

template<class T>
class const_n_pointer<T,1,C_TYPE,true>
{
        const T* p_p;
        const size_t* p_st;
        const tree_vec* p_v;
public:
        const_n_pointer(const T* p, const size_t* st, const tree_vec* v) : p_p(p), p_st(st), p_v(v) {}
        const T& operator[] (const size_t i) const
        {
                if( i >= p_v->n )
                        OUT_OF_RANGE( "const_n_pointer::operator[]" );
                return *(p_p + i);
        }
};

//
//! multi_arr: generic class for allocating multidimensional arrays.
//! A typical example of its use could be:
//!
//!     multi_arr<double,3> arr; // define a placeholder for the array
//!                       // the first argument is the type of data it holds
//!                       // the second argument is the number of dimensions
//!                       // (between 2 and 6)
//!
//!     arr.alloc(3,4,2); // this will allocate a 3x4x2 block of doubles
//!                       // memory will be allocated by new[] so will
//!                       // be initialized by the constructor, if present
//! 
//!     multi_arr<double,3> arr(3,4,2); // shorthand for the above
//!
//! The following is an alternative way of allocating the array. It is
//! very similar to the pre-multi_arr way of allocating arrays. In ARPA_TYPE
//! arrays this will help you save memory since only the data elements
//! that are really needed will be allocated. In C_TYPE allocation, the
//! smallest rectangular block will be allocated that can hold all the data.
//! This will use more memory in return for somewhat improved CPU speed.
//! Tests carried out in 2007 showed that the speed advantage of C_TYPE arrays
//! was only 1 to 2 percent. Hence the memory savings were deemed more important
//! and ARPA_TYPE arrays were made the default. However, C_TYPE arrays are
//! guaranteed to be compatible with C code, so these should be used if they
//! are meant to be passed on to C library routines. The example below allocates
//! a triangular matrix.
//!
//!     arr.reserve(3);
//!     for( int i=0, i < 3; ++i )
//!     {
//!       arr.reserve( i, i+1 ); // note that size does not need to be constant!
//!       for( int j=0, j < i+1; ++j )
//!         arr.reserve( i, j, j+1 );
//!     }
//!     arr.alloc();
//!
//! these are plausible ways to use the multi_arr class:
//!        
//!     arr.invalidate(); // this will set float or double arrays to all SNaN
//!                       // it will set any other type array to all 0xff bytes.
//!     arr.zero();       // this will set the array to all zero
//!
//!     arr[0][0][0] = 1.;
//!     arr[0][0][1] = 2.;
//!     double x = arr[0][0][0]*arr[0][0][1];
//!
//!     arr.state_do(io,lgRead); // this will write/read the array to/from a state file
//!                       // depending on the boolean flag lgRead; io should point to a
//!                       // file that is already opened for binary writing/reading
//!                       // the file will remain open after access. this allows you to
//!                       // use state_do() for several arrays in a row
//!
//!     multi_arr<double,2,C_TYPE> a(10,10); // allocate C_TYPE array
//!     C_library_routine( a.data(), ... );  // and call C library routine with it
//!
//!     arr.clear();      // this will deallocate the array
//!                       // the destructor will also automatically deallocate
//!
//! the multi_arr class comes with iterators that allow you to speed up memory access
//! even further. using iterators as shown below will generally speed up the code
//! significantly since it avoids calculating the array index over and over inside the
//! body of the loop. especially in tight loops over arrays with high dimension this can
//! become a significant overhead! a const_iterator is also supplied for read-only access,
//! but no reverse_iterators. you can define and initialize an iterator as follows
//!
//!     multi_arr<double,3>::iterator p = arr.begin(n,k);
//!
//! the notation multi_arr