// © Copyright 2008 Vladimir Frolov, Moscow State University Graphics & Media Lab
#pragma once

#include <valarray>
#include <vector>
#include <queue>

#include "MGML_ERROR.h"

using namespace MGML_ERROR;

namespace MGML_MEMORY
{


// заменитель std::vector
template<class T>
class vector
{
public:

	void	 push_back(const T& val);
	const T& pop_back();
	void	 erase(int k);
  void	 reserve(int n);
  void	 resize(int n);

	typedef T*		 iterator;
	typedef const T* const_iterator;

	inline const T& operator[](int i) const { return elems[i];}
	inline       T& operator[](int i){ return elems[i];}
	inline const const_iterator begin() const {return elems;}
	inline const const_iterator end()   const {return elems + N;}
	
	inline unsigned int size() const {return N;}

	vector();
	vector(int n);
	~vector();

	vector& operator=(const vector&);

protected:
	enum {mem_stride      = 16,
		  unrolling_param = 4};

	void elemcpy(T* to,T* from,int N);


	T* elems;
	unsigned int N;
	unsigned int mem_tos;

};

// простой вектор, занимает мешьше памяти
template<class T>
class simple_vector
{
public:
	typedef T*		 iterator;
	typedef const T* const_iterator;

	simple_vector();
	simple_vector(int n);
	~simple_vector();

	//void	 push_back(const T& val);
	//void	 reserve(int n);
	void		 resize(int n);   // обрезаем лишние элементы.
	unsigned int size() const 
	{
		return N;
	}

	inline const T operator[](int i) const { return elems[i];}
	inline       T& operator[](int i){ return elems[i];}
	inline const const_iterator begin() const {return elems;}
	inline const const_iterator end()   const {return elems + N;}
	inline const T* get_pointer() const {return elems;}

	simple_vector& operator=(const simple_vector&);

protected:
	enum {unrolling_param = 4};

	void elemcpy(T* to,const T* from,int N);


	T* elems;
	unsigned int N;
};



template<int n,class lword = UCHAR>
class Bit_Array
{
public:

 Bit_Array(){}
 Bit_Array(size_t n_): arr(n_*n){}

 class ref
 {
	 int i1;
	 Bit_Array<n,lword> *v;
 public:
	 ref(){}
	 ref(Bit_Array<n,lword>* v1,int j1): v(v1),i1(j1) {}

	 inline ref& operator=(lword val)
		{ v->set(i1,val); return *this; }

	 inline operator lword() const {return v->get(i1);}
 };

 inline ref operator[](int i)
  { return ref(this,i); }

 inline const ref operator[](int i) const
  { return ref(const_cast< Bit_Array<n,lword>* >(this),i); }

 //inline const lword operator[](int i) const
 // { return get(i); }

 inline const size_t size() const 
  {return (arr.size() >> 1);}

 inline void push_back(lword arg)
 {
	for(int i=0;i<n;i++) 
		arr.push_back((arg >> (n-i-1)) & 0x1);
 }

 inline lword pop_back()
 {
	 register lword res = get((arr.size() >> 1) - 1);
	 for(int i=0;i<n;i++) 
		arr.pop_back();
	 return res;
 }

 inline const lword get(int i) const
 {
	register lword res = arr[i*n] << (n-1);
	for(int j=1;j<n;j++)
		  res |= arr[(i*n)+j] << (n-j-1);
	return res;
 }

 inline void set(int index,lword arg)
 {
	for(int i=0;i<n;i++) 
		arr[n*index + i] = (arg >> (n-i-1)) & 0x1;
 }

protected:
	std::vector<bool> arr;
};



template<class lword>
class Bit_Array<2,lword>
{
public:
	Bit_Array()
		{if(sizeof(lword) != 1) throw 0;}
	Bit_Array(int _n): arr(_n/4 + 1)
		{if(sizeof(lword) != 1) throw 0; _size = _n;}
 
 class ref
 {
	 int i1;
	 Bit_Array<2,lword> *v;
 public:
	 ref(){}
	 ref(Bit_Array<2,lword>* v1,int j1): v(v1),i1(j1) {}

	 inline ref& operator=(lword val)
		{ v->set(i1,val); return *this; }

	 inline operator lword() const {return v->get(i1);}
 };

 inline const size_t size() const 
  {return (_size);}

 inline ref operator[](int i)
  { return ref(this,i); }

 inline const ref operator[](int i) const
  { return ref(const_cast< Bit_Array<2,lword>* >(this),i); }

 void resize(int _n) // не отлажено, может быть источником ошибок
 {
	 int n2 = _n/4 + 1;
	 int min_ = MIN<int>(arr.size(),n2);
     arr.resize(min_);
	 _size = _n;
 }
 /*inline void push_back(lword arg)
 {
	 register size_t vsize = arr.size() << 2;
	 register int rem = vsize - size();
	 if(rem == 0)
		 arr.push_back(arg);
	 else
	 {
		 switch(rem)
		 {
			case 1: arr.[vsize].e3 = arg; break;
			case 2: arr.[vsize].e2 = arg; break;
			case 3: arr.[vsize].e1 = arg; break;
		 };
	 }
 }*/

 inline const lword get(int i) const
 {
	register int ind = (i >> 2);
	register int rem = i - (ind << 2);
	register int res; 
	switch(rem)
	{
		case 0:res = arr[ind].e1; break;
		case 1:res = arr[ind].e2; break;
		case 2:res = arr[ind].e3; break;
		case 3:res = arr[ind].e4; break;
	};
	return res;
 }

 inline void set(int index,lword arg)
 {
	register int ind = (index >> 2);
	register int rem = index - (ind << 2);
	switch(rem)
	{
		case 0: arr[ind].e1 = arg; break;
		case 1: arr[ind].e2 = arg; break;
		case 2: arr[ind].e3 = arg; break;
		case 3: arr[ind].e4 = arg; break;
	};
 }

private:

	struct st_elem
	{
		union
		{
			lword  num;
			struct
			{
				lword e1 : 2;
				lword e2 : 2;
				lword e3 : 2;
				lword e4 : 2;
			};
		};
	};

	std::valarray<st_elem> arr;
	size_t _size;

};





////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////  Matrix  ///////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////

template
<
	int n,						// размерность матрицы (2D 3D 4D итд.)
	class T,					// тип элементов
	class TR,					// тип возвр. ссылок (выбирать по умолчанию T&)
	template <class> class vec	// класс, использующийся для представления вектора
>
class Matrix
{

public:

	template<int m> class slice;
	template<int m> class slice_iter;
	template<int m> class cslice_iter;

	Matrix(const size_t _s[n]) 
	{
		for(int i=0,len = 1;i<n;i++)
			s[i] = _s[i];
		v = new vec<T>(len_arr());
		set_strides();
	}
	Matrix(const Matrix& m2) 
	{
		for(int i=0;i<n;i++)  s[i] = 0;
		v = new vec<T>(m2.const_array().size());
		*this = m2;
		set_strides();
	}
	Matrix()
	{
		for(int i=0,len = 1;i<n;i++) 
			s[i] =0;
		v = new vec<T>(0);
		set_strides();
	}
	~Matrix(){ delete v;}

	Matrix<n,T,TR,vec>& operator=(const Matrix<n,T,TR,vec>& m2);

	inline size_t len_arr() const
	{
		size_t len = 1;
		for(int i=0;i<n;i++) len*=s[i];
		return len;
	}

	
	void resize(size_t N[n]) 
	{ 
		delete v;
		size_t len = 1;
		for(int i=0;i<n;i++)
		{
			s[i] = N[i];
			len*=s[i];
		}
		v = new vec<T>(len); 
		set_strides();
	}

	template<int coord> 
	inline size_t size() const {return s[coord];}
	
	template<int coord>
	inline size_t stride() const {return strides[coord];}

	size_t offset(const size_t _s[n]) const
	{
		register int offset = 0;
		for(int i=0;i<n-1;i++) 
			offset += strides[i]*_s[i];
		offset += _s[n-1];

		return offset;
	}

	inline vec<T>& array() {return *v;}
	inline const vec<T>& const_array() const {return *v;}

	inline slice_iter<n-1> row(size_t i)
	 { return slice_iter<n-1>(v,slice<n-1>(i*stride<0>(),s+1));  }

	inline cslice_iter<n-1> row(size_t i) const
	 { return cslice_iter<n-1>(v,slice<n-1>(i*stride<0>(),s+1));  }

	inline TR operator()(const size_t _s[n]) { return array()[offset(_s)]; }
	inline T  operator()(const size_t _s[n]) const { return const_array()[offset(_s)]; }

	inline slice_iter<n-1> operator[](size_t i) {return row(i);}
	inline cslice_iter<n-1> operator[](size_t i) const {return row(i);}

	template<int m>
	class slice
	{
	public:

		inline slice(size_t start_,const size_t* ss): _start(start_)
		 { for(int i=0;i<m;i++) sz[i]=ss[i]; }

		inline size_t stride() const
		{
			size_t res=1;
			for(int i=1;i<m;i++)
				res *= sz[i];
			return res;
		}

		inline size_t size() const
		{
			size_t res=1;
			for(int i=0;i<m;i++) res *= sz[i];
			return res;
		}

		inline size_t start() const { return _start; }
		inline const size_t* get_sz() const {return sz;}

	protected:
		size_t _start;
		size_t sz[m];

	};

	template<int m>
	class slice_iter
	{
		typedef slice_iter<m-1> TypeRef;

	public:
		inline slice_iter(vec<T>* _v,const slice<m>& _s): v(_v),s(_s) {}

		inline slice_iter<m>& operator++(){curr++;return *this;}
		inline TypeRef operator[](size_t i) {return ref(curr = i);}
		inline TypeRef operator()(size_t i) {return ref(curr = i);}
		inline TypeRef operator*(){return ref(curr);}
		inline size_t size() {return s.size();}
	
	protected:
		vec<T>*  v;
		slice<m> s;
		size_t	 curr;

		inline TypeRef ref(size_t i) const 
		 { return TypeRef(v,slice<m-1>(i*s.stride() + s.start(),s.get_sz()+1) ); }
	};

	template<>
	class slice_iter<1>
	{

	public:
		slice_iter(vec<T>* _v,const slice<1>& _s): v(_v),s(_s) {}
		
		inline slice_iter<1>& operator++(){curr++;return *this;}
		inline TR operator[](size_t i) {return ref(curr = i);}
		inline TR operator()(size_t i) {return ref(curr = i);}
		inline TR operator*(){return ref(curr);}
		inline size_t size() {return s.size();}

	protected:
		vec<T>*  v;
		slice<1> s;
		size_t	 curr;

		inline TR ref(size_t i) const { return (*v)[s.start() + i*s.stride()]; }
	};

	template<int m>
	class cslice_iter
	{
		typedef cslice_iter<m-1> TypeRef;

	public:
		inline cslice_iter(const vec<T>* _v,const slice<m>& _s): v(_v),s(_s) {}

		inline cslice_iter<m>& operator++(){curr++;return *this;}
		inline TypeRef operator[](size_t i) {return ref(curr = i);}
		inline TypeRef operator()(size_t i) {return ref(curr = i);}
		inline TypeRef operator*(){return ref(curr);}
		inline size_t size() {return s.size();}
	
	protected:
		const vec<T>*  v;
		slice<m> s;
		size_t	 curr;

		inline TypeRef ref(size_t i) const 
		 { return TypeRef(v,slice<m-1>(i*s.stride() + s.start(),s.get_sz()+1) ); }
	};

	template<>
	class cslice_iter<1>
	{

	public:
		cslice_iter(const vec<T>* _v,const slice<1>& _s): v(_v),s(_s) {}
		
		inline slice_iter<1>& operator++(){curr++;return *this;}
		inline T operator[](size_t i) {return ref(curr = i);}
		inline T operator()(size_t i) {return ref(curr = i);}
		inline T operator*(){return ref(curr);}
		inline size_t size() {return s.size();}

	protected:
		const vec<T>*  v;
		slice<1> s;
		size_t	 curr;

		inline T ref(size_t i) const { return (*v)[s.start() + i*s.stride()]; }
	};

protected:

	inline void set_strides()
	{
		for(int i=0;i<n-1;i++)
		{
			strides[i] = 1;
			for(int j=i+1;j<n;j++) 
				strides[i]*= s[j];
		}
	}

	vec<T>*		 v; 
	size_t		 s[n];
	size_t       strides[n-1];

};




template<class T>
class FastList
{
public:

  class Elem
  {
  public:

    friend class FastList;
  protected:
    T m_data;
    Elem* m_next;
  };

  class iterator
  {
  public:

    iterator() {m_elem = 0;}
    iterator(Elem* a_elem) 
    {
      m_elem = a_elem;
    }

    bool operator==(const iterator& a_it)
    { return (a_it.m_elem == m_elem); }

    bool operator!=(const iterator& a_it)
    { return !(*this == a_it); }

    iterator operator++()
    {
      ASSERT(m_elem!=0);
      m_elem = m_elem->m_next;
      return iterator(m_elem);
    }

    iterator operator++(int)
    {
      iterator prev = *this;
      ++(*this);
      return prev;
    }

    T* operator->() const 
    {
      ASSERT(m_elem!=0);
      return &(m_elem->m_data);
    }

    T& operator*() 
    { 
      ASSERT(m_elem!=0);
      return m_elem->m_data;
    }


  protected:
    Elem* m_elem;
  };


  class const_iterator
  {
  public:

    const_iterator() {m_elem = 0;}
    const_iterator(const Elem* a_elem) 
    {
      m_elem = a_elem;
    }

    bool operator==(const const_iterator& a_it)
    { return (a_it.m_elem == m_elem); }

    bool operator!=(const const_iterator& a_it)
    { return !(*this == a_it); }

    iterator operator++()
    {
      ASSERT(m_elem!=0);
      m_elem = m_elem->m_next;
      return const_iterator(m_elem);
    }

    iterator operator++(int)
    {
      const_iterator prev = *this;
      ++(*this);
      return prev;
    }

    T* operator->() const 
    {
      ASSERT(m_elem!=0);
      return &(m_elem->m_data);
    }

    T& operator*() 
    { 
      ASSERT(m_elem!=0);
      return m_elem->m_data;
    }


  protected:
    const Elem* m_elem;
  };


  FastList() { m_head = 0; }

  inline bool Empty() {return (m_head == 0);}
  inline void Reset() { m_head = 0;} // you must Reset all your lists after Free call, or delete them 

  inline iterator begin() { return iterator(m_head);}
  inline iterator end() { return iterator(0);}

  inline const_iterator begin() const { return iterator(m_head);}
  inline const_iterator end() const { return iterator(0);}


  inline void push_front(const T& a_data)
  {
    Elem *newElem = m_freeList;
    ASSERT(newElem!=0);
    m_freeList = newElem->m_next;
    newElem->m_next = m_head;
    newElem->m_data = a_data;
    m_head = newElem;
  }

  inline void pop_front()
  {
    if(m_head == 0)
      return;

    Elem* next = m_head->m_next;
    m_head->m_next = m_freeList;
    m_freeList = m_head;
    m_head = next;   
  }

  void clear()
  {
    while(m_head != 0)
      pop_front();
  }

  size_t size() const // а нехрен меня вызывать
  {
    size_t count = 0;
    Elem *newElem = m_head;
    while( newElem )
    {
      ++count;
      newElem = newElem->m_next;
    };
    return count;                 
  }

  static void Allocate(int a_size)
  {
    m_memory = new Elem[a_size];

    for(int i=0;i<a_size-1;i++)
      m_memory[i].m_next = m_memory + i + 1;
    m_memory[a_size-1].m_next = 0;

    m_freeList = m_memory;
  }

  static void Free() 
  { 
    delete [] m_memory;  
    m_memory = 0;
  }

protected:
  static Elem* m_memory;
  static Elem* m_freeList;

  Elem* m_head;
};

template<class T> typename FastList<T>::Elem* FastList<T>::m_memory = 0;
template<class T> typename FastList<T>::Elem* FastList<T>::m_freeList = 0;





};




////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////  Binary_Tree  ///////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////


template<class NodeData, class Allocator = std::allocator<NodeData> >
class Binary_Tree
{
public:


	struct Node
	{
		Node()
		{
			m_Data = 0xFFFFFFFF;
			data = NodeData();
		}
		

		void SetLeftOffset(unsigned int a_Left)
		{
			register int mask = a_Left & 0x3fffffff;
			register int m_Data2 = m_Data & 0xc0000000;
			m_Data = m_Data2 | mask;
		}


		unsigned int GetLeftOffset() const 
			{ return m_Data & 0x3FFFFFFF;}

		unsigned int GetRightOffset() const 
			{ return GetLeftOffset()+1;}


		bool HasLeftChild() const
			{ return !((m_Data & 0x80000000) >> 31); } // one bit for right sub-tree empty flag

		bool HasRightChild() const
			{ return !((m_Data & 0x40000000) >> 30); } // one bit for right sub-tree empty flag

		bool Leaf() const 
			{ return !HasLeftChild() & !HasRightChild();}  // one bit to know - is there a leaf?


		void SetLeftChildEmpty(bool arg)
		{
			register int mask = ((int)arg) << 31;
			register int m_Data2 = m_Data & 0x7fffffff;
			m_Data = m_Data2 | mask;
		}

		void SetRightChildEmpty(bool arg)
		{
			register int mask = ((int)arg) << 30;
			register int m_Data2 = m_Data & 0xbfffffff;
			m_Data = m_Data2 | mask;
		}


		NodeData data;
	protected:
		int m_Data;
		
	};
	

	Binary_Tree()
	{
		nodes_array.resize(0);
		nodes_top = 1;
		//Node* pThisNode = getCurrNode();
	}

	size_t GetRootOffset() {return 0;}

	const Node& GetNodeByOffset(size_t offset) const
	{
		ASSERT(offset < nodes_top);
		return nodes_array[offset];
	}

	Node& GetNodeByOffset(size_t offset)
	{
		ASSERT(offset < nodes_top);
		return nodes_array[offset];
	}

	void Print(std::ostream& out)
	{
		out << "Tree: " << endl;
		print(out,GetRootOffset(),0);
	}

	struct compare_simple
	{
		int operator()(const NodeData& n1, const NodeData& n2)
			{ return n1 < n2;}
	};

	template<class Compare>
	void Build(NodeData* a_DataArray, unsigned int a_Num, Compare cmp)
	{
		nodes_array.resize(a_Num*3-1);
		
		tmp_arr.resize(0);
		tmp_arr.reserve(a_Num);

		for(size_t i=0;i<a_Num;i++)
			tmp_arr.push_back(a_DataArray[i]);

		std::sort(tmp_arr.begin(),tmp_arr.end(),cmp);

//		std::vector<NodeData>::const_iterator p;
//		for(p=tmp_arr.begin();p!=tmp_arr.end();++p)
//		{
//			std::cout << *p <<" ";
//		}
//		std::cout << std::endl;

		subdivide(0,tmp_arr.size()-1,GetRootOffset());

		Reorganize_memory_layout_BFS(nodes_queue);

		tmp_arr.resize(0);
	}

	std::vector<Node>& GetNodesArray() { return nodes_array; }

protected:
	std::vector<Node> nodes_array;
	size_t nodes_top;

	std::vector<NodeData> tmp_arr;
	
	// for advanced memory layout
	std::queue<size_t> nodes_queue;

	Binary_Tree(const Binary_Tree& rhs){}
	Binary_Tree& operator=(const Binary_Tree& rhs)
		{return *this;}

	void print(std::ostream& out,size_t offs,int deep)
	{
		for(int i=0;i<deep;i++)
			out << " ";

		const Node& node = GetNodeByOffset(offs);
		out <<"offs: " << offs << ", data: " << node.data << std::endl;

		if(node.HasLeftChild())
			print(out,node.GetLeftOffset(),deep+2);

		if(node.HasRightChild())
			print(out,node.GetRightOffset(),deep+2);
		
	}

	size_t NewNodePair()
	{
		if(nodes_top+2 >= nodes_array.size())
			throw Out_Of_Range("Binary_Tree, too many elements");

		Node nodes[2];
		nodes_array[nodes_top+0] = nodes[0];
		nodes_array[nodes_top+1] = nodes[1];
		nodes_top+=2;

		return nodes_top-2;
	}

	void subdivide(size_t begin, size_t end, size_t curr_offset)
	{
		size_t size = end-begin+1;
		if(size == 1)
		{
			Node* pThisNode = &GetNodeByOffset(curr_offset);
			pThisNode->data = tmp_arr[begin];
			pThisNode->SetRightChildEmpty(true);
			pThisNode->SetLeftChildEmpty(true);

			return;
		}
		else if(size == 2)
		{
			Node* pThisNode    = &GetNodeByOffset(curr_offset);
			size_t left_offset = NewNodePair();
			pThisNode->data    = tmp_arr[begin+1];
			pThisNode->SetRightChildEmpty(true);
			pThisNode->SetLeftChildEmpty(false);
			pThisNode->SetLeftOffset((unsigned int)left_offset);

			Node* pLeftNode	   = &GetNodeByOffset(left_offset);
			pLeftNode->data    = tmp_arr[begin];
			pLeftNode->SetRightChildEmpty(true);
			pLeftNode->SetLeftChildEmpty(true);

			return;
		}
		//	throw Error("Error in Binary_Tree building algorithm");
	
		size_t	median		= begin + (size-1)/2;
		size_t	begin_left	= begin;
		size_t	end_left	= median-1;
		size_t	begin_right	= median+1;
		size_t	end_right	= end;

		Node* pThisNode     = &GetNodeByOffset(curr_offset);
		size_t left_offset  = NewNodePair();
		size_t right_offset = left_offset+1;

		
		pThisNode->SetRightChildEmpty(false);
		pThisNode->SetLeftChildEmpty(false);
		pThisNode->SetLeftOffset((unsigned int)left_offset);
		pThisNode->data = tmp_arr[median];

		subdivide(begin_left, end_left,	left_offset);
		subdivide(begin_right,end_right,right_offset);
		
	}

	void BFS_traversal(size_t a_Offset, std::queue<size_t>& a_Queue)
	{
		std::queue<size_t> s_Queue;

		s_Queue.push(a_Offset);
		while(!s_Queue.empty())
		{

			const Node& node = GetNodeByOffset(s_Queue.front());
			a_Queue.push(s_Queue.front());
			s_Queue.pop();


			if(node.HasLeftChild())
			{
				s_Queue.push(node.GetLeftOffset());
				s_Queue.push(node.GetRightOffset());
			}
		}
	}

	void Reorganize_memory_layout_BFS(std::queue<size_t>& a_Queue)
	{
		size_t offs = GetRootOffset();
		BFS_traversal(offs,a_Queue);
		std::vector<Node>	nodes_array2(nodes_top);
		std::vector<size_t>	new_by_old_offs(nodes_top);

		ASSERT(a_Queue.size() == nodes_top);

		size_t i = 0;
		while(!a_Queue.empty())
		{
			size_t new_offs = i;
			size_t old_offs = a_Queue.front();
			a_Queue.pop();

			new_by_old_offs[old_offs] = (unsigned int)new_offs;
			i++;
		}

		for(size_t i = 0;i<nodes_top;i++)
		{
			size_t new_offs = new_by_old_offs[i]; 
			nodes_array2[new_offs] = nodes_array[i];

			size_t old_left_offs = nodes_array[i].GetLeftOffset();
			size_t new_left_offs;
			if(nodes_array[i].HasLeftChild())
				new_left_offs = new_by_old_offs[old_left_offs];
			else
				new_left_offs = 0xFFFFFFFF;
			
			nodes_array2[new_offs].SetLeftOffset((unsigned int)new_left_offs);
		}

		nodes_array.resize(nodes_array2.size());
		for(size_t i=0;i<nodes_array2.size();i++)
			nodes_array[i] = nodes_array2[i];
	}
};









////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////  vector  ///////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////

template<class T>
MGML_MEMORY::vector<T>::vector()
{
 elems = 0;
 N  = 0;
 mem_tos =0;
}

template<class T>
MGML_MEMORY::vector<T>::vector(int n)
{
 elems = new T [n];
 N  = n;
 mem_tos =n;
}


template<class T>
MGML_MEMORY::vector<T>::~vector()
{
 delete [] elems;
}

template<class T>
MGML_MEMORY::vector<T>& MGML_MEMORY::vector<T>::operator=(const vector<T>& rhs)
{
	if(size() < rhs.size())
		reserve(rhs.size());

	int n1 = size();

	N = rhs.size();

	for(UINT i=0;i<N;i++)
		elems[i] = rhs[i];
	
	return *this;
}

template<class T>
void MGML_MEMORY::vector<T>::reserve(int n)
{
	if(elems)
	{
		if((unsigned)n < N)
			throw Error("strange requirement to reserve less memory than already reserved"); // вообще это скорее всего ошибка

		T* temp = new T[size()];
		elemcpy(temp,elems,size());
		delete [] elems;
		elems = new T [n];
		elemcpy(elems,temp,size());
		delete [] temp; // strange errors
	}
	else 
		elems = new T [n];

    mem_tos = n;
}


template<class T>
void MGML_MEMORY::vector<T>::resize(int n)
{
	if(elems)
	{
		T* temp = new T[size()];
		elemcpy(temp,elems,size());
		delete [] elems;
		elems = new T [n];
		elemcpy(elems,temp,size());
		delete [] temp; // strange errors
	}
	else 
		elems = new T [n];

  N = n;
  mem_tos = n;
}


template<class T>
void MGML_MEMORY::vector<T>::push_back(const T& val)
{
	if(size() + 1 >= mem_tos)
      reserve(mem_stride+mem_tos);

	elems[N] = val;
	N++;
}


template<class T>
const T& MGML_MEMORY::vector<T>::pop_back()
{
	N--;
	return elems[N];
}

template<class T>
void MGML_MEMORY::vector<T>::erase(int k)
{
	if(elems)
	{
		elemcpy(elems+k,elems+k+1,N-k);
		N--;
	}
	else
		throw 0;
}


template<class T>
void MGML_MEMORY::vector<T>::elemcpy(T* to,T* from,int N)
{
 int N_remain = N%unrolling_param;
 int N_main_iter = (N -N_remain);
 
 for(int i=0;i<N_main_iter;i+=unrolling_param)
	 for(int j=0;j<unrolling_param;j++)
		to[i+j] = from[i+j];

 for(int i=N_main_iter;i<N_main_iter+N_remain;i++)
	 to[i] = from[i];

}

////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////  simple_vector  ///////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////
template<class T>
MGML_MEMORY::simple_vector<T>::simple_vector(int n)
{
 elems = new T [n];
 N  = n;
}

template<class T>
MGML_MEMORY::simple_vector<T>::simple_vector()
{
 elems = 0;
 N     = 0;
}

template<class T>
MGML_MEMORY::simple_vector<T>::~simple_vector()
{
	delete [] elems;
}


template<class T>
void MGML_MEMORY::simple_vector<T>::elemcpy(T* to,const T* from,int N)
{
 int N_remain = N%unrolling_param;
 int N_main_iter = N - N_remain;
 
 for(int i=0;i<N_main_iter;i+=unrolling_param)
	 for(int j=0;j<unrolling_param;j++)
		to[i+j] = from[i+j];

 for(int i=N_main_iter;i<N_main_iter+N_remain;i++)
	 to[i] = from[i];

}


template<class T>
void MGML_MEMORY::simple_vector<T>::resize(int n)
{
	if(elems)
	{
		
		if((unsigned)n == N) return; 
		N = n;
		T* temp = new T[n];
		elemcpy(temp,elems,n);
		delete [] elems;
		elems = new T [n];
		elemcpy(elems,temp,n);
		delete [] temp;
	}
	else 
	{
		elems = new T [n];
		N = n;
	}
}

template<class T>
MGML_MEMORY::simple_vector<T>& MGML_MEMORY::simple_vector<T>::operator=(const simple_vector& rhs)
{
	if(size() != rhs.size())
		resize(rhs.size());

	elemcpy(elems,rhs.get_pointer(),size());
	return *this;
}
/*
template<class T>
void MGML_MEMORY::simple_vector<T>::reserve(int n)
{
	if(elems)
	{
		if((unsigned)n < N) return; // вообще это скорее всего ошибка

		T* temp = new T[size()];
		elemcpy(temp,elems,size());
		delete [] elems;
		elems = new T [n];
		elemcpy(elems,temp,size());
		delete [] temp;
	}
	else 
		elems = new T [n];

    mem_tos = n;
}

template<class T>
void MGML_MEMORY::simple_vector<T>::push_back(const T& val)
{
	if(size() + 1 >= size())
      reserve(1+size());

	elems[N] = val;
	N++;
}*/

////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////  Matrix  ///////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////
//======================================================================================================
// 
//======================================================================================================
template
<
	int n,
	class T,					// тип элементов
	class TR,					// тип возвр. елементов (T или T&)
	template <class> class vec	// класс, использующийся для представления вектора
>
MGML_MEMORY::Matrix<n,T,TR,vec>& MGML_MEMORY::Matrix<n,T,TR,vec>::operator=(const Matrix<n,T,TR,vec>& m2)
{
	if(v->size() < m2.const_array().size())	 
	{
		delete v;
		int n1 = m2.const_array().size();
		v = new vec<T>(m2.const_array().size());
	}
	
	*v  = m2.const_array();

	for(int i=0;i<n;i++)
		s[i] = m2.s[i];

	set_strides();

	return *this;
}