container.hpp

#ifndef _CONTAINER_HPP
#define _CONTAINER_HPP


#include <stdlib.h>
#include <stddef.h>
#include <string.h>
#include <memory>
#include <new>
#include <iterator>
#include <cmath>
#include <assert.h>
#include <limits>
#include "inmost_common.h"

//#define DEBUGMEM

#define NEW_CHUNKS
#define __VOLATILE volatile
//#define OUT_OF_RANGE

#if defined(_OPENMP)
#include <omp.h>
#endif

//TODO
// 1. change to uniform size_type instead of size_t, make it INMOST_DATA_ENUM_TYPE
/*
template<class element, class T1> struct isInputRandomIterators
{
	static void constraints(T1 a, T1 b) { ++a; (void)a++; a==a; a!=a; a-b; }
	isInputRandomIterators() { void(*p)(T1,T1) = constraints; (void)p; }
};

template<class element, class T1> struct isInputForwardIterators
{
	static void constraints(T1 a) {++a; (void)a++; }
	isInputForwardIterators() { void(*p)(T1) = constraints; (void)p; }
};
*/

namespace INMOST
{
	//add more templates as needed
	template<class T> struct make_unsigned;
	template<> struct make_unsigned<char> {typedef unsigned char type;};
	template<> struct make_unsigned<short> {typedef unsigned short type;};
	template<> struct make_unsigned<int> {typedef unsigned int type;};
	template<> struct make_unsigned<long> {typedef unsigned long type;};
	template<> struct make_unsigned<long long> {typedef unsigned long long type;};
	template<> struct make_unsigned<unsigned char> {typedef unsigned char type;};
	template<> struct make_unsigned<unsigned short> {typedef unsigned short type;};
	template<> struct make_unsigned<unsigned int> {typedef unsigned int type;};
	template<> struct make_unsigned<unsigned long> {typedef unsigned long type;};
	template<> struct make_unsigned<unsigned long long> {typedef unsigned long long type;};

//#define USE_OPTIMIZED_ARRAY_ALLOCATION
#if defined(PACK_ARRAY)
#pragma pack(push,r1,4)
#endif
	// notice: array class would not properly support classes that contain self-references
	//         like std::map
	// notice: next class shell have to implement same algorithms as array
	template<typename element>//, typename enumerator = unsigned int>
	class array
	{
	public:
		typedef unsigned size_type;
		typedef make_unsigned<size_type>::type uenum;
		template<typename etype>
		class _iterator
		{
		private:
			etype * e;
		public:
			typedef etype * pointer;
			typedef etype & reference;
			typedef etype value_type;
			typedef ptrdiff_t difference_type;
			typedef std::random_access_iterator_tag iterator_category;
			_iterator():e(NULL){}
			_iterator(etype * i):e(i){}
			_iterator(const _iterator & other){e = other.e;}
			~_iterator() {};
			_iterator operator -(size_t n) const { return _iterator(e-n); }
			_iterator & operator -=(size_t n) { e-=n; return *this; }
			_iterator operator +(size_t n) const { return _iterator(e+n); }
			_iterator & operator +=(size_t n) { e+=n; return *this; }
			_iterator & operator ++(){ ++e; return *this;}
			_iterator operator ++(int){ return _iterator(e++); }
			_iterator & operator --(){ --e; return *this; }
			_iterator operator --(int){ return _iterator(e--); }
			ptrdiff_t operator -(const _iterator & other) const {return e-other.e;}
			etype & operator *() { return *e; }
			const etype & operator *() const { return *e; }
			etype * operator ->() { return e; }
			_iterator & operator =(_iterator const & other) { e = other.e; return *this; }
			bool operator ==(const _iterator & other) const { return e == other.e;}
			bool operator !=(const _iterator & other) const { return e != other.e;}
			bool operator <(const _iterator & other) const { return e < other.e;}
			bool operator >(const _iterator & other) const { return e > other.e;}
			bool operator <=(const _iterator & other) const { return e <= other.e;}
			bool operator >=(const _iterator & other) const { return e >= other.e;}
			operator void *() {return static_cast<void *> (e);}
			operator const void *() const {return const_cast<const void *> (e);}
		};
		typedef _iterator<element> iterator;
		typedef _iterator<const element> const_iterator;
		template<typename etype>
		class _reverse_iterator
		{
		private:
			etype * e;
		public:
			typedef etype * pointer;
			typedef etype & reference;
			typedef etype value_type;
			typedef ptrdiff_t difference_type;
			typedef std::random_access_iterator_tag iterator_category;
			_reverse_iterator():e(NULL){}
			_reverse_iterator(etype * i):e(i){}
			_reverse_iterator(const _reverse_iterator & other){e = other.e;}
			~_reverse_iterator() {};
			_reverse_iterator operator -(size_t n) const { return _reverse_iterator(e+n); }
			_reverse_iterator & operator -=(size_t n) { e+=n; return *this; }
			_reverse_iterator operator +(size_t n) const {return _reverse_iterator(e-n); }
			_reverse_iterator & operator +=(size_t n) { e-=n; return *this; }
			_reverse_iterator & operator ++(){ --e; return *this;}
			_reverse_iterator operator ++(int){ return _reverse_iterator(e--); }
			_reverse_iterator & operator --(){ ++e; return *this; }
			_reverse_iterator operator --(int){ return _reverse_iterator(e++); }
			ptrdiff_t operator -(const _reverse_iterator & other) const {return other.e-e;}
			etype & operator *() { return *e; }
			const etype & operator *() const { return *e; }
			etype * operator ->() { return e; }
			_reverse_iterator & operator =(_reverse_iterator const & other) { e = other.e; return *this;}
			bool operator ==(const _reverse_iterator & other) const { return e == other.e;}
			bool operator !=(const _reverse_iterator & other) const { return e != other.e;}
			bool operator <(const _reverse_iterator & other) const { return e < other.e;}
			bool operator >(const _reverse_iterator & other) const { return e > other.e;}
			bool operator <=(const _reverse_iterator & other) const { return e <= other.e;}
			bool operator >=(const _reverse_iterator & other) const { return e >= other.e;}
			operator void *() {return static_cast<void *> (e);}
			operator const void *() const {return static_cast<const void *> (e);}
		};
		typedef _reverse_iterator<element> reverse_iterator;
		typedef _reverse_iterator<const element> const_reverse_iterator;
	private:

		element * m_arr;
		size_type m_size;

		//notice: push_back, pop_back, insert, erase are optimized for current growth_formula
		__INLINE static size_type growth_formula(size_type future_size)
		{
			uenum v = static_cast<uenum>(future_size);
			v--;
			v|= (v>>1);
			v|= (v>>2);
			v|= (v>>4);
			v|= (v>>8);
			v|= (v>>16);
			//TODO produces compiler warning
			//if( sizeof(size_type) > 4 ) v|= (v>>32); //must be compile time brench
			v++;
			return static_cast<size_type>(v);
		}
	public:
		__INLINE element * data() {return m_arr;}
		__INLINE const element * data() const {return m_arr;}
		array() {m_arr = NULL; m_size = 0;}
		array(size_type n,element c = element())
		{
			m_size = n;
			void * tmp = malloc(sizeof(element)*growth_formula(m_size));
#if defined(DEBUGMEM)
			if( tmp == NULL ) {printf("%s:%d malloc return NULL\n",__FILE__,__LINE__);}
#endif
			m_arr = static_cast<element *>(tmp);
			assert(m_arr != NULL);
			for(size_type i = 0; i < m_size; i++) new (m_arr+i) element(c);
		}
		template<class InputIterator>
		array(InputIterator first, InputIterator last)
		{
			//isInputForwardIterators<element,InputIterator>();
			m_size = static_cast<size_type>(std::distance(first,last));
			void * tmp = malloc(sizeof(element)*growth_formula(m_size));
#if defined(DEBUGMEM)
			if( tmp == NULL ) {printf("%s:%d malloc return NULL\n",__FILE__,__LINE__);}
#endif
			m_arr = static_cast<element *>(tmp);
			assert(m_arr != NULL);
			{
				size_type i = 0;
				InputIterator it = first;
				while(it != last) new (m_arr+i++) element(*it++);
			}
		}
		array(const array & other)
		{
			m_size = other.m_size;
			if( m_size )
			{
				void * tmp = malloc(sizeof(element)*growth_formula(m_size));
#if defined(DEBUGMEM)
				if( tmp == NULL ) {printf("%s:%d malloc return NULL\n",__FILE__,__LINE__);}
#endif
				m_arr = static_cast<element *>(tmp);
				assert(m_arr != NULL);
			}
			else m_arr = NULL;
			for(size_type i = 0; i < m_size; i++) new (m_arr+i) element(other.m_arr[i]);
		}
		~array()
		{
			for(size_type i = 0; i < m_size; i++) m_arr[i].~element();
			if( m_arr != NULL ) free(m_arr);
			m_arr = NULL;
			m_size = 0;
		}
		__INLINE const element & operator [] (size_type n) const
		{
			assert(n < m_size);
			return m_arr[n];
		}
		__INLINE element & operator [] (size_type n)
		{
			assert(n < m_size);
			return m_arr[n];
		}
		__INLINE const element & at (size_type n) const
		{
			assert(n < m_size);
			return m_arr[n];
		}
		__INLINE element & at (size_type n)
		{
			assert(n < m_size);
			return m_arr[n];
		}
		__INLINE element & at_safe (size_type n)
		{
			if( n >= m_size ) resize(n+1);
			return m_arr[n];
		}
		array & operator =(array const & other)
		{
			if( this != &other )
			{
				for(size_type i = 0; i < m_size; i++) m_arr[i].~element();
				if(m_arr != NULL )
				{
					free(m_arr);
					m_arr = NULL;
					m_size = 0;
				}
				if( other.m_arr != NULL )
				{
					m_size = other.m_size;
					void * tmp = malloc(sizeof(element)*growth_formula(m_size));
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d malloc return NULL\n",__FILE__,__LINE__);}
#endif
					m_arr = static_cast<element *>(tmp);
					assert(m_arr != NULL);
					if( m_arr ) memcpy(m_arr,other.m_arr,sizeof(element)*m_size);
				}
			}
			return *this;
		}
		void push_back(const element & e)
		{
#if !defined(USE_OPTIMIZED_ARRAY_ALLOCATION)
			//unoptimized variant
			if( m_size+1 > growth_formula(m_size) )
			{
				void * tmp = realloc(m_arr,sizeof(element)*growth_formula(++m_size));
#if defined(DEBUGMEM)
				if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif				
				m_arr = static_cast<element *>(tmp);
			}
			else ++m_size;
#else
			//optimized for current growth_formula
			if( m_size < 2 )
			{
				void * tmp = realloc(m_arr,sizeof(element)*(++m_size));
#if defined(DEBUGMEM)
				if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif		
				m_arr = static_cast<element *>(tmp);
			}
			else if( ((m_size+1) & (m_size-1)) == 1 )
			{
				void * tmp = realloc(m_arr,sizeof(element)*(m_size++ << 1));
#if defined(DEBUGMEM)
				if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif		
				m_arr = static_cast<element *>(tmp);
			}
			else m_size++;
#endif
			assert(m_arr != NULL);
			new (m_arr+m_size-1) element(e);
		}
		void pop_back()
		{
			assert(m_arr != NULL);
			m_arr[m_size--].~element();
			if( m_size > 0)
			{
#if !defined(USE_OPTIMIZED_ARRAY_ALLOCATION)
				//unoptimized variant
				size_type gf = growth_formula(m_size);
				if( m_size+1 > gf )
				{
					void * tmp = realloc(m_arr,sizeof(element)*gf);
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif		
					m_arr = static_cast<element *>(tmp);
				}
#else
				if( ((m_size+1) & (m_size-1)) == 1 || m_size == 1)
				{
					void * tmp = realloc(m_arr,sizeof(element)*m_size);
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif		
					m_arr = static_cast<element *>(tmp);
				}
#endif
				assert(m_arr != NULL);
			}
			else
			{
				free(m_arr);
				m_arr = NULL;
			}
		}
		__INLINE element & back()
		{
			assert(m_arr != NULL);
			assert(m_size > 0);
			return m_arr[m_size-1];
		}
		__INLINE const element & back() const
		{
			assert(m_arr != NULL);
			assert(m_size > 0);
			return m_arr[m_size-1];
		}
		__INLINE element & front()
		{
			assert(m_arr != NULL);
			assert(m_size > 0);
			return m_arr[0];
		}
		__INLINE const element & front() const
		{
			assert(m_arr != NULL);
			assert(m_size > 0);
			return m_arr[0];
		}
		__INLINE size_type capacity() { return m_size; }
		__INLINE bool empty() const { if( m_size ) return false; return true; }
		void resize(size_type n, element c = element() )
		{
			size_type oldsize = m_size;
			m_size = n;
			for(size_type i = m_size; i < oldsize; i++) m_arr[i].~element(); //delete elements, located over the size
			if( m_size > 0 )
			{
				if( growth_formula(oldsize) != growth_formula(m_size) )
				{
					void * tmp = realloc(m_arr,sizeof(element)*growth_formula(m_size));
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					m_arr = static_cast<element *>(tmp);
				}
				assert(m_arr != NULL);
				for(size_type i = oldsize; i < m_size; i++) new (m_arr+i) element(c); //initialize extra entities
			}
			else
			{
				free(m_arr);
				m_arr = NULL;
			}
		}
		__INLINE size_type size() const {return m_size;}
		__INLINE size_type capacity() const {return growth_formula(m_size);}
		void clear()
		{
			for(size_type i = 0; i < m_size; i++) m_arr[i].~element();
			m_size = 0;
			if( m_arr ) free(m_arr);
			m_arr = NULL;
		}
		void swap(array<element> & other)
		{
			element * t_m_arr = m_arr;
			size_type t_m_size = m_size;
			m_arr = other.m_arr;
			m_size = other.m_size;
			other.m_arr = t_m_arr;
			other.m_size = t_m_size;
		}
		__INLINE iterator begin() { return m_arr; }
		__INLINE iterator end() { return m_arr+m_size; }
		__INLINE const_iterator begin() const { return m_arr; }
		__INLINE const_iterator end() const { return m_arr+m_size; }
		__INLINE reverse_iterator rbegin() { return reverse_iterator(m_arr+m_size-1); }
		__INLINE reverse_iterator rend() { return reverse_iterator(m_arr-1); }
		__INLINE const_reverse_iterator rbegin() const { return const_reverse_iterator(m_arr+m_size-1); }
		__INLINE const_reverse_iterator rend() const { return const_reverse_iterator(m_arr-1); }
		iterator erase(iterator pos)
		{
			ptrdiff_t d = pos-begin();
			ptrdiff_t s = iterator(m_arr+m_size-1)-pos;
			(*pos).~element();
			m_size--;
			if( m_size > 0 )
			{
				if( s > 0 ) memmove(m_arr+d,m_arr+d+1,sizeof(element)*s);
				//unoptimized variant
#if !defined(USE_OPTIMIZED_ARRAY_ALLOCATION)
				size_type gf = growth_formula(m_size);
				if( m_size+1 > gf )
				{
					void * tmp = realloc(m_arr,sizeof(element)*gf);
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					m_arr = static_cast<element *>(tmp);
				}
#else
				if( ((m_size+1) & (m_size-1)) == 1 || m_size == 1)
				{
					void * tmp = realloc(m_arr,sizeof(element)*m_size);
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					m_arr = static_cast<element *>(tmp);
				}
#endif
				assert(m_arr != NULL);
			}
			else
			{
				free(m_arr);
				m_arr = NULL;
			}
			return m_arr+d;
		}
		iterator erase(iterator b, iterator e)
		{
			ptrdiff_t d = b-begin();
			ptrdiff_t s = end()-e;
			ptrdiff_t n = e-b;
			for(iterator i = b; i != e; i++) (*i).~element();
			m_size -= n;
			if( m_size > 0 )
			{
				if( s > 0 ) memmove(m_arr+d,m_arr+d+1,sizeof(element)*s);
				size_type gf = growth_formula(m_size);
				if( m_size+n > gf )
				{
					void * tmp = realloc(m_arr,sizeof(element)*gf);
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					m_arr = static_cast<element *>(tmp);
				}
				assert(m_arr != NULL);
			}
			else
			{
				free(m_arr);
				m_arr = NULL;
			}
			return m_arr+d;
		}
		iterator insert(iterator pos, const element & x)
		{
			if( static_cast<void *>(pos) == NULL)
			{
				assert(m_arr == NULL);
				pos = iterator(m_arr = static_cast<element *>(malloc(sizeof(element))));
#if defined(DEBUGMEM)
				if( m_arr == NULL ) {printf("%s:%d malloc return NULL\n",__FILE__,__LINE__);}
#endif
				assert(m_arr != NULL);
			}
			ptrdiff_t d = pos-begin();
			ptrdiff_t s = end()-pos;

			//unoptimized variant
#if !defined(USE_OPTIMIZED_ARRAY_ALLOCATION)
			if( m_size+1 > growth_formula(m_size) )
			{
				void * tmp = realloc(m_arr,sizeof(element)*growth_formula(++m_size));
#if defined(DEBUGMEM)
				if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
				m_arr = static_cast<element *>(tmp);
			}
			else ++m_size;
#else
			//optimized for current growth_formula
			if( m_size < 2 )
			{
				void * tmp = realloc(m_arr,sizeof(element)*(++m_size));
#if defined(DEBUGMEM)
				if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
				m_arr = static_cast<element *>(tmp);
			}
			else if( ((m_size+1) & (m_size-1)) == 1 )
			{
				void * tmp = realloc(m_arr,sizeof(element)*(m_size++ << 1));
#if defined(DEBUGMEM)
				if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
				m_arr = static_cast<element *>(tmp);
			}
			else ++m_size;
#endif
			assert(m_arr != NULL);
			if( s > 0 ) memmove(m_arr+d+1,m_arr+d,sizeof(element)*s);
			new (m_arr+d) element(x);
			return m_arr+d;
		}
		void insert(iterator pos, size_type n, const element & x)
		{
			if( n > 0 )
			{
				if( static_cast<void *>(pos) == NULL)
				{
					assert(m_arr == NULL);
					pos = iterator(m_arr = static_cast<element *>(malloc(sizeof(element))));
#if defined(DEBUGMEM)
					if( m_arr == NULL ) {printf("%s:%d malloc return NULL\n",__FILE__,__LINE__);}
#endif
					assert(m_arr != NULL);
				}
				ptrdiff_t d = pos-begin();
				ptrdiff_t s = end()-pos;

				if( m_size+n > growth_formula(m_size) )
				{
					void * tmp = realloc(m_arr,sizeof(element)*growth_formula(m_size+n));
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					m_arr = static_cast<element *>(tmp);
				}
				m_size+=n;

				assert(m_arr != NULL);
				if( s > 0 ) memmove(m_arr+d+n,m_arr+d,sizeof(element)*s);
				for(size_type i = 0; i < n; i++) new (m_arr+d+i) element(x);
			}
		}
		template <class InputIterator>
		void insert(iterator pos, InputIterator first, InputIterator last)
		{
			size_type n = static_cast<size_type>(std::distance(first,last));
			if( n > 0 )
			{
				if( static_cast<void *>(pos) == NULL)
				{
					assert(m_arr == NULL);
					pos = iterator(m_arr = static_cast<element *>(malloc(sizeof(element))));
#if defined(DEBUGMEM)
					if( m_arr == NULL ) {printf("%s:%d malloc return NULL\n",__FILE__,__LINE__);}
#endif
					assert(m_arr != NULL);
				}
				ptrdiff_t d = pos-begin();
				ptrdiff_t s = end()-pos;

				if( m_size+n > growth_formula(m_size) )
				{
					void * tmp = realloc(m_arr,sizeof(element)*growth_formula(m_size+n));
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					m_arr = static_cast<element *>(tmp);
				}
				m_size+=n;

				assert(m_arr != NULL);
				if( s > 0 ) memmove(m_arr+d+n,m_arr+d,sizeof(element)*s);
				{
					InputIterator it = first;
					size_type i = 0;
					while(it != last) new (m_arr+d+i++) element(*it++);
				}
			}
		}
		template <class InputIterator>
		void replace(iterator m_first, iterator m_last, InputIterator first, InputIterator last)
		{
			assert( m_size >= 0 );
			ptrdiff_t n = static_cast<ptrdiff_t>(std::distance(first,last));
			if( static_cast<void *>(m_first) == NULL && m_arr == NULL)
			{
				assert(m_arr == NULL);
				m_first = m_last = iterator(m_arr = static_cast<element *>(malloc(sizeof(element))));
#if defined(DEBUGMEM)
				if( m_arr == NULL ) {printf("%s:%d malloc return NULL\n",__FILE__,__LINE__);}
#endif
				assert(m_arr != NULL);
			}
			ptrdiff_t q = m_last-m_first;
			ptrdiff_t d = m_first-iterator(m_arr);
			ptrdiff_t s = iterator(m_arr+m_size)-m_last;
			for(iterator it = m_first; it != m_last; it++) (*it).~element();
			if( n-q != 0 )
			{
				size_type gf = growth_formula(m_size+static_cast<size_type>(n-q));
				if( gf != growth_formula(m_size) )
				{
					void * tmp = realloc(m_arr,sizeof(element)*gf);
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					m_arr = static_cast<element *>(tmp);
				}
				m_size+=static_cast<size_type>(n-q);
				if( s > 0 ) memmove(m_arr+d+n,m_arr+d+q,sizeof(element)*s);
			}
			{
				InputIterator it = first;
				size_type i = 0;
				while(it != last) new (m_arr+d+i++) element(*it++);
			}
		}
		template <class InputIterator>
		void assign(InputIterator first, InputIterator last)
		{
			replace(begin(),end(),first,last);
		}
		template<class> friend class shell;
	};
#if defined(PACK_ARRAY)
#pragma pack(pop,r1)
#endif
	template<typename element>
	class shell
	{
	public:
		typedef typename array<element>::size_type size_type;
		template<typename dtype>
		class _iterator
		{
		private:
			dtype * e;
		public:
			typedef dtype * pointer;
			typedef dtype & reference;
			typedef dtype value_type;
			typedef ptrdiff_t difference_type;
			typedef std::random_access_iterator_tag iterator_category;
			_iterator():e(NULL){}
			_iterator(dtype * i):e(i){}
			_iterator(const _iterator & other){e = other.e;}
			~_iterator() {};
			_iterator operator -(size_t n) const { return _iterator(e-n); }
			_iterator & operator -=(size_t n) { e-=n; return *this; }
			_iterator operator +(size_t n) const { return _iterator(e+n); }
			_iterator & operator +=(size_t n) { e+=n; return *this; }
			_iterator & operator ++(){ ++e; return *this;}
			_iterator operator ++(int){ return _iterator(e++); }
			_iterator & operator --(){ --e; return *this; }
			_iterator operator --(int){ return _iterator(e--); }
			ptrdiff_t operator -(const _iterator & other) const {return e-other.e;}
			dtype & operator *() { return *e; }
			const dtype & operator *() const { return *e; }
			dtype * operator ->() { return e; }
			_iterator & operator =(_iterator const & other) { e = other.e; return *this; }
			bool operator ==(const _iterator & other) { return e == other.e;}
			bool operator !=(const _iterator & other) { return e != other.e;}
			bool operator <(const _iterator & other) { return e < other.e;}
			bool operator >(const _iterator & other) { return e > other.e;}
			bool operator <=(const _iterator & other) { return e <= other.e;}
			bool operator >=(const _iterator & other) { return e >= other.e;}
			operator void *() {return static_cast<void *> (e);}
		};
		typedef _iterator<element> iterator;
		typedef _iterator<const element> const_iterator;
		template<typename dtype>
		class _reverse_iterator
		{
		private:
			dtype * e;
		public:
			typedef dtype * pointer;
			typedef dtype & reference;
			typedef dtype value_type;
			typedef ptrdiff_t difference_type;
			typedef std::random_access_iterator_tag iterator_category;
			_reverse_iterator():e(NULL){}
			_reverse_iterator(dtype * i):e(i){}
			_reverse_iterator(const _reverse_iterator & other){e = other.e;}
			~_reverse_iterator() {};
			_reverse_iterator operator -(size_t n) const { return _reverse_iterator(e+n); }
			_reverse_iterator & operator -=(size_t n) { e+=n; return *this; }
			_reverse_iterator operator +(size_t n) const {return _reverse_iterator(e-n); }
			_reverse_iterator & operator +=(size_t n) { e-=n; return *this; }
			_reverse_iterator & operator ++(){ --e; return *this;}
			_reverse_iterator operator ++(int){ return _reverse_iterator(e--); }
			_reverse_iterator & operator --(){ ++e; return *this; }
			_reverse_iterator operator --(int){ return _reverse_iterator(e++); }
			ptrdiff_t operator -(const _reverse_iterator & other) const {return other.e-e;}
			dtype & operator *() { return *e; }
			const dtype & operator *() const { return *e; }
			dtype * operator ->() { return e; }
			_reverse_iterator & operator =(_reverse_iterator const & other) { e = other.e; return *this;}
			bool operator ==(const _reverse_iterator & other) { return e == other.e;}
			bool operator !=(const _reverse_iterator & other) { return e != other.e;}
			bool operator <(const _reverse_iterator & other) { return e < other.e;}
			bool operator >(const _reverse_iterator & other) { return e > other.e;}
			bool operator <=(const _reverse_iterator & other) { return e <= other.e;}
			bool operator >=(const _reverse_iterator & other) { return e >= other.e;}
			operator void *() {return static_cast<void *> (e);}
		};
		typedef _reverse_iterator<element> reverse_iterator;
		typedef _reverse_iterator<const element> const_reverse_iterator;
	private:
		element ** m_arr;
		size_type * m_size;
		element * local_link;
		size_type local_size;
		bool fixed;
	public:
		__INLINE element * data() {return *m_arr;}
		__INLINE const element * data() const {return *m_arr;}
		shell() :local_link(NULL), local_size(0), m_arr(NULL), m_size(NULL), fixed(false) { }
		shell(array<element> & arr) //dynamic
			:local_link(NULL), local_size(0), m_arr(&arr.m_arr), m_size(&arr.m_size), fixed(false) {}
		shell(element * link, size_type size) //fixed
			:local_link(NULL),local_size(0), m_arr(&local_link), m_size(&local_size), fixed(true)
		{
			*m_arr = link;
			*m_size = size;
		}
		shell(const shell & other)
			:local_link(NULL), local_size(0), m_arr(&local_link), m_size(&local_size),fixed(other.fixed)
		{
			if( fixed )
			{
				*m_size = *other.m_size;
				*m_arr = *other.m_arr;
			}
			else
			{
				m_size = other.m_size;
				m_arr = other.m_arr;
			}
		}
		~shell()
		{
		}
		__INLINE const element & operator [] (size_type n) const
		{
			assert(n < *m_size );
			return *((*m_arr)+n);
		}
		__INLINE element & operator [] (size_type n)
		{
			assert(n < *m_size );
			return *((*m_arr)+n);
		}
		__INLINE const element & at (size_type n) const
		{
			assert(n < *m_size );
			return *((*m_arr)+n);
		}
		__INLINE element & at (size_type n)
		{
			assert(n < *m_size );
			return *((*m_arr)+n);
		}
		shell & operator =(shell const & other)
		{
			fixed = other.fixed;
			if( fixed )
			{
				m_size = &local_size;
				m_arr = &local_link;
				*m_size = *other.m_size;
				*m_arr = *other.m_arr;
			}
			else
			{
				m_size = other.m_size;
				m_arr = other.m_arr;
			}
			return *this;
		}
		void push_back(const element & e)
		{
			assert( !fixed ); // array size is fixed
#if !defined(USE_OPTIMIZED_ARRAY_ALLOCATION)
			//unoptimized variant
			if( (*m_size)+1 > array<element>::growth_formula(*m_size) )
			{
				void * tmp = realloc(*m_arr,sizeof(element)*array<element>::growth_formula(++(*m_size)));
#if defined(DEBUGMEM)
				if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif				
				*m_arr = static_cast<element *>(tmp);
			}
			else ++(*m_size);
#else
			//optimized for current growth_formula
			if( *m_size < 2 )
			{
				void * tmp = realloc(*m_arr,sizeof(element)*(++(*m_size)));
#if defined(DEBUGMEM)
				if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
				*m_arr = static_cast<element *>(tmp);
			}
			else if( (((*m_size)+1) & ((*m_size)-1)) == 1 )
			{
				void * tmp = realloc(*m_arr,sizeof(element)*((*m_size)++ << 1));
#if defined(DEBUGMEM)
				if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
				*m_arr = static_cast<element *>(tmp);
			}
			else ++(*m_size);
#endif
			assert((*m_arr) != NULL);
			new ((*m_arr)+(*m_size)-1) element(e);
		}
		void pop_back()
		{
			assert( !fixed ); // array size is fixed
			assert((*m_arr) != NULL);
			(*m_arr)[(*m_size)--].~element();
			if( (*m_size) > 0 )
			{
#if !defined(USE_OPTIMIZED_ARRAY_ALLOCATION)
				//unoptimized variant
				size_type gf = array<element>::growth_formula(*m_size);
				if( (*m_size)+1 > gf )
				{
					void * tmp = realloc(m_arr,sizeof(element)*gf);
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					*m_arr = static_cast<element *>(tmp);
				}
#else
				if( (((*m_size)+1) & ((*m_size)-1)) == 1 || (*m_size) == 1)
				{
					void * tmp = realloc(*m_arr,sizeof(element)*(*m_size));
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					*m_arr = static_cast<element *>(tmp);
				}
#endif
				assert( (*m_arr) != NULL );
			}
			else
			{
				free(*m_arr);
				(*m_arr) = NULL;
			}
		}
		__INLINE element & back()
		{
			assert(*m_arr != NULL);
			assert(*m_size > 0 );
			return (*m_arr)[(*m_size)-1];
		}
		__INLINE const element & back() const
		{
			assert(*m_arr != NULL);
			assert(*m_size > 0 );
			return (*m_arr)[(*m_size)-1];
		}
		__INLINE element & front()
		{
			assert(*m_arr != NULL);
			assert(*m_size > 0 );
			return (*m_arr)[0];
		}
		__INLINE const element & front() const
		{
			assert(*m_arr != NULL);
			assert(*m_size > 0 );
			return (*m_arr)[0];
		}
		__INLINE size_type capacity() { return array<element>::growth_formula(*m_size); }
		__INLINE bool empty() const { if( *m_size ) return false; return true; }
		void resize(size_type n, element c = element() )
		{
			assert( !fixed ); // array size is fixed
			size_type oldsize = *m_size;
			*m_size = n;
			for(size_type i = *m_size; i < oldsize; i++) (*m_arr)[i].~element(); //delete elements, located over the size
			if( *m_size > 0 )
			{
				if( array<element>::growth_formula(oldsize) != array<element>::growth_formula(*m_size) )
				{
					void * tmp = realloc(*m_arr,sizeof(element)*array<element>::growth_formula(*m_size));
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					*m_arr = static_cast<element *>(tmp);
				}
				assert( (*m_arr) != NULL );
				for(size_type i = oldsize; i < *m_size; i++) new ((*m_arr)+i) element(c); //initialize extra entities
			}
			else
			{
				free(*m_arr);
				*m_arr = NULL;
			}
		}
		__INLINE size_type size() const {return *m_size;}
		void clear()
		{
			assert( !fixed ); // array size is fixed
			for(size_type i = 0; i < *m_size; i++) (*m_arr)[i].~element();
			*m_size = 0;
			if( *m_arr ) free(*m_arr);
			*m_arr = NULL;
		}
		void swap(shell<element> & other)
		{
			element * t_m_arr = *m_arr;
			size_type t_m_size = *m_size;
			*m_arr = *other.m_arr;
			*m_size = *other.m_size;
			*other.m_arr = t_m_arr;
			*other.m_size = t_m_size;
		}
		__INLINE iterator begin() { return *m_arr; }
		__INLINE iterator end() { return *m_arr+(*m_size); }
		__INLINE const_iterator begin() const { return *m_arr; }
		__INLINE const_iterator end() const { return *m_arr+(*m_size); }
		__INLINE reverse_iterator rbegin() { return reverse_iterator(*m_arr+(*m_size)-1); }
		__INLINE reverse_iterator rend() { return reverse_iterator(*m_arr-1); }
		__INLINE const_reverse_iterator rbegin() const { return const_reverse_iterator(*m_arr+(*m_size)-1); }
		__INLINE const_reverse_iterator rend() const { return const_reverse_iterator(*m_arr-1); }
		iterator erase(iterator pos)
		{
			assert( !fixed ); // array size is fixed
			ptrdiff_t d = pos-begin();
			ptrdiff_t s = iterator(*m_arr+(*m_size)-1)-pos;
			(*pos).~element();
			(*m_size)--;
			if( (*m_size) > 0 )
			{
				if( s > 0 )memmove(*m_arr+d,*m_arr+d+1,sizeof(element)*s);
#if !defined(USE_OPTIMIZED_ARRAY_ALLOCATION)
				size_type gf = array<element>::growth_formula(*m_size);
				if( (*m_size)+1 > gf )
				{
					void * tmp = realloc(*m_arr,sizeof(element)*gf);
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					*m_arr = static_cast<element *>(tmp);
				}
#else
				if( (((*m_size)+1) & ((*m_size)-1)) == 1 || (*m_size) == 1)
				{
					void * tmp = realloc(*m_arr,sizeof(element)*(*m_size));
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					*m_arr = static_cast<element *>(tmp);
				}
#endif
				assert((*m_arr) != NULL);
			}
			else
			{
				free(*m_arr);
				*m_arr = NULL;
			}
			return (*m_arr)+d;
		}
		iterator erase(iterator b, iterator e)
		{
			assert( !fixed ); // array size is fixed
			ptrdiff_t d = b-begin();
			ptrdiff_t s = end()-e;
			ptrdiff_t n = e-b;
			for(iterator i = b; i != e; i++) (*i).~element();
			(*m_size) -= n;
			if( (*m_size) > 0 )
			{
				if( s > 0 ) memmove(*m_arr+d,*m_arr+d+n,sizeof(element)*s);
				size_type gf = array<element>::growth_formula(*m_size);
				if( (*m_size)+n > gf )
				{
					void * tmp = realloc(*m_arr,sizeof(element)*gf);
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					*m_arr = static_cast<element *>(tmp);
				}
				assert((*m_arr) != NULL);
			}
			else
			{
				free(*m_arr);
				*m_arr = NULL;
			}
			return (*m_arr)+d;
		}
		iterator insert(iterator pos, const element & x)
		{
			assert( !fixed ); // array size is fixed
			if( static_cast<void *>(pos) == NULL )
			{
				assert((*m_arr) == NULL);
				pos = iterator((*m_arr) = static_cast<element *>(malloc(sizeof(element))));
#if defined(DEBUGMEM)
				if( *m_arr == NULL ) {printf("%s:%d malloc return NULL\n",__FILE__,__LINE__);}
#endif
				assert((*m_arr) != NULL);
			}
			ptrdiff_t d = pos-begin();
			ptrdiff_t s = end()-pos;

#if !defined(USE_OPTIMIZED_ARRAY_ALLOCATION)
			//unoptimized variant
			if( (*m_size)+1 > array<element>::growth_formula(*m_size) )
			{
				void * tmp = realloc(*m_arr,sizeof(element)*array<element>::growth_formula(++(*m_size)));
#if defined(DEBUGMEM)
				if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
				*m_arr = static_cast<element *>(tmp);
			}
			else ++(*m_size);
#else
			//optimized for current growth_formula
			if( *m_size < 2 )
			{
				void * tmp = realloc(*m_arr,sizeof(element)*(++(*m_size)));
#if defined(DEBUGMEM)
				if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
				*m_arr = static_cast<element *>(tmp);
			}
			else if( (((*m_size)+1) & ((*m_size)-1)) == 1 )
			{
				void * tmp = realloc(*m_arr,sizeof(element)*((*m_size)++ << 1));
#if defined(DEBUGMEM)
				if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
				*m_arr = static_cast<element *>(tmp);
			}
			else ++(*m_size);
#endif
			assert((*m_arr) != NULL);
			if( s > 0 ) memmove((*m_arr)+d+1,(*m_arr)+d,sizeof(element)*s);
			new ((*m_arr)+d) element(x);
			return (*m_arr)+d;
		}
		void insert(iterator pos, size_type n, const element & x)
		{
			assert( !fixed ); // array size is fixed
			if( n > 0 )
			{
				if( static_cast<void *>(pos) == NULL)
				{
					assert((*m_arr) == NULL);
					pos = iterator((*m_arr) = static_cast<element *>(malloc(sizeof(element))));
#if defined(DEBUGMEM)
					if( *m_arr == NULL ) {printf("%s:%d malloc return NULL\n",__FILE__,__LINE__);}
#endif
					assert((*m_arr) != NULL);
				}
				ptrdiff_t d = pos-iterator(*m_arr);
				ptrdiff_t s = iterator((*m_arr)+(*m_size))-pos;

				if( (*m_size)+n > array<element>::growth_formula(*m_size) )
				{
					void * tmp = realloc(*m_arr,sizeof(element)*array<element>::growth_formula((*m_size)+n));
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					*m_arr = static_cast<element *>(tmp);
				}
				(*m_size)+=n;

				assert((*m_arr) != NULL);
				if( s > 0 ) memmove((*m_arr)+d+n,(*m_arr)+d,sizeof(element)*s);
				for(size_type i = 0; i < n; i++) new ((*m_arr)+d+i) element(x);
			}
		}
		template <class InputIterator>
		void insert(iterator pos, InputIterator first, InputIterator last)
		{
			assert( !fixed ); // array size is fixed
			ptrdiff_t n = static_cast<ptrdiff_t>(std::distance(first,last));
			if( n > 0 )
			{
				if( static_cast<void *>(pos) == NULL)
				{
					assert((*m_arr) == NULL);
					pos = iterator((*m_arr) = static_cast<element *>(malloc(sizeof(element))));
#if defined(DEBUGMEM)
					if( *m_arr == NULL ) {printf("%s:%d malloc return NULL\n",__FILE__,__LINE__);}
#endif
					assert((*m_arr) != NULL);
				}
				ptrdiff_t d = pos-iterator(*m_arr);
				ptrdiff_t s = iterator((*m_arr)+(*m_size))-pos;


				if( (*m_size)+n > array<element>::growth_formula(*m_size) )
				{
					void * tmp = realloc(*m_arr,sizeof(element)*array<element>::growth_formula((*m_size)+static_cast<size_type>(n)));
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					*m_arr = static_cast<element *>(tmp);
				}
				(*m_size)+=static_cast<size_type>(n);

				assert((*m_arr) != NULL);
				if( s > 0 ) memmove((*m_arr)+d+n,(*m_arr)+d,sizeof(element)*s);
				{
					InputIterator it = first;
					size_type i = 0;
					while(it != last) new ((*m_arr)+d+i++) element(*it++);
				}
			}
		}
		template <class InputIterator>
		void replace(iterator m_first, iterator m_last, InputIterator first, InputIterator last)
		{
			ptrdiff_t n = static_cast<ptrdiff_t>(std::distance(first,last));
			if( static_cast<void *>(m_first) == NULL)
			{
				assert((*m_arr)==NULL);
				m_first = m_last = iterator((*m_arr) = static_cast<element *>(malloc(sizeof(element))));
#if defined(DEBUGMEM)
				if( *m_arr == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
				assert((*m_arr)!=NULL);
			}
			ptrdiff_t q = m_last-m_first;
			ptrdiff_t d = m_first-iterator(*m_arr);
			ptrdiff_t s = iterator((*m_arr)+(*m_size))-m_last;
			for(iterator it = m_first; it != m_last; it++) (*it).~element();
			if( n-q != 0 )
			{
				assert( !fixed ); // array size is fixed
				size_type gf = array<element>::growth_formula((*m_size)+static_cast<size_type>(n-q));
				if( gf != array<element>::growth_formula(*m_size) )
				{
					void * tmp = realloc(*m_arr,sizeof(element)*gf);
#if defined(DEBUGMEM)
					if( tmp == NULL ) {printf("%s:%d realloc return NULL\n",__FILE__,__LINE__);}
#endif
					*m_arr = static_cast<element *>(tmp);
				}
				(*m_size)+=static_cast<size_type>(n-q);
			}
			if( s > 0 )
				memmove((*m_arr)+d+n,(*m_arr)+d+q,sizeof(element)*s);
			{
				InputIterator it = first;
				size_type i = 0;
				while(it != last) new ((*m_arr)+d+i++) element(*it++);
			}
		}
		template <class InputIterator>
		void assign(InputIterator first, InputIterator last)
		{
			replace(begin(),end(),first,last);
		}
	};



	template <typename IndType,typename ValType>
	class sparse_data
	{
	public:
		typedef unsigned enumerate;
		static const enumerate prealloc = 4;
		typedef struct pair_t
		{
			IndType first;
			ValType second;
			pair_t() :first(0),second(0.0) {}
			pair_t(IndType first, ValType second) :first(0), second(0.0) {}
		} pair;
		typedef pair * iterator;
		typedef const iterator const_iterator;
	private:
		static size_t comparator(const void * pa, const void *