#include "binMatrix.h"



//--------------------------------------------------------------------------
// -- this part should be modified if the type packet is modified
//    3,7,8,255 are magic number related to the number of bits per packet
const packetSize = 8;
const fullP = 255;
const packet p2 [8]= {   1,   2,   4,   8,  16,  32,  64, 128 };
const packet p2c[8]= { 254, 253, 251, 247, 239, 223, 191, 127 }; // 255-p2
const packet p2l[8]= {   1,   3,   7,  15,  31,  63, 127, 255 }; // sum p2
const packet p2r[8]= { 254, 252, 248, 240, 224, 192, 128,   0 }; // 255-p2l
//--------------------------------------------------------------------------




////////////////////////
// Array Constructors //
////////////////////////

/*inline*/ binArray::binArray(int l, int m, int n, int jOff, packet* d/**=0**/) :
	ArrayP(l, m, ceilDB(n+jOff), d), offJ(jOff), realN(n) { 
}


/*inline*/ binArray::binArray(int m, int n, int jOff, packet* d/**=0**/) :
	ArrayP(ceilDB(n+jOff), m, ceilDB(n+jOff), d), offJ(jOff), realN(n) { 
}


/*inline*/ binArray::binArray(int m, int n, packet* d/**=0**/) :
	ArrayP(ceilDB(n), m, ceilDB(n), d), offJ(0), realN(n) { 
}


/*inline*/ binArray::binArray(int n, packet* d/**=0**/) :
	ArrayP(ceilDB(n), 1, ceilDB(n), d), offJ(0), realN(n) { 
}


/*inline*/ binArray::binArray(const binArray& a) :
	ArrayP(a), offJ(a.offJ), realN(a.realN) { 
}


/*inline*/ binArray::binArray() : ArrayP(), offJ(0), realN(0) { 
}


/*inline*/ binArray::~binArray() { 
}



//////////////////////
// Arrays selectors //
//////////////////////


/*inline*//*??*/ int binArray::m() const
{	return ArrayP::m(); 
}


int binArray::n() const
{	return realN; 
}


// subArray (without space reallocation)

/*inline*/ binArray	binArray::s(int i, int m, int j, int n)	const
{
#ifdef _CHECK_INDEX
  if ( i<0 )
    throw whereOutOfRange("binArray binArray::s(int,int,int,int) const",
			  "i",i,0,this->m()-1);
  if ( m<0 || i+m>this->m() )
    throw whereOutOfRange("binArray binArray::s(int,int,int,int) const",
			  "i+m",i+m,0,this->m()-1);
  if ( j<0 )
    throw whereOutOfRange("binArray binArray::s(int,int,int,int) const",
			  "j",j,0,this->n()-1);
  if ( n<0 || j+n>this->n() )
    throw whereOutOfRange("binArray binArray::s(int,int,int,int) const",
			  "j+n",j+n,0,this->n()-1);
#endif

  return binArray(l(), m, n, moduloB(j+offJ), x()+i*l()+floorDB(j+offJ));
}


/*inline*/ binArray	binArray::s(int i/**=0**/, int m/**=1**/, int j/**=0**/) const
{
	#ifdef _CHECK_INDEX
	if ( i<0 )
		throw whereOutOfRange(
		  "binArray binArray::s(int=0,int=1,int=0) const","i",i,0,this->m()-1);
	if ( m<0 || i+m>this->m() )
		throw whereOutOfRange(
			"binArray binArray::s(int=0,int=1,int=0) const",
			"i+m",i+m,0,this->m()-1);
	if ( j<0 || j>=n() )
		throw whereOutOfRange(
		  "binArray binArray::s(int=0,int=1,int=0) const","j",j,0,this->n()-1);
	#endif
	return binArray(l(), m, n()-j, moduloB(j+offJ), x()+i*l()+floorDB(j+offJ));
}


// subMatrix (with space reallocation)

/*inline*/ binMatrix binArray::_(int i, int m, int j, int n) const
{	return binMatrix(s(i, m, j, n)); 
}

/*inline*/ binMatrix binArray::_(int i/**=0**/, int m/**=1**/, int j/**=0**/) const
{	return binMatrix(s(i, m, j, n()-j)); 
}


// subscripting

/*inline*/ boolean binArray::operator () (int i, int j) const
{
	#ifdef _CHECK_INDEX
	if ( i<0 || i>=m() )
		throw whereOutOfRange(
			"boolean binArray::operator()(int,int) const","i",i,0,m()-1);
	if ( j<0 || j>=n() )
		throw whereOutOfRange(
			"boolean binArray::operator()(int,int) const","j",j,0,n()-1);
	#endif
	return ( ( ArrayP::operator()(i,floorDB(j+offJ)) & p2[moduloB(j+offJ)] ) ?
			 true : false );
}


/*inline*/ boolean binArray::operator () (int j) const
{	return binArray::operator()(0,j); 
}


/*inline*/ char binArray::block (int i, int j) const
{
	#ifdef _CHECK_INDEX
	if ( i<0 || i>=m() )
		throw whereOutOfRange(
			"char binArray::block(int,int) const","i",i,0,m()-1);
	if ( j<0 || j>=ArrayP::n() )
		throw whereOutOfRange(
			"boolean binArray::operator()(int,int) const",
			"j",j,0,ArrayP::n()-1);
	#endif
	if ( (j && j < ArrayP::n()-1) || 
		 (!j && !offJ && realN >= packetSize) ||
		 (j==ArrayP::n()-1 && !moduloB(realN+offJ) && realN >= packetSize) )
		return ( ArrayP::operator()(i,j) );
	else if (!j )
		if ( offJ )
			if ( realN+offJ >= packetSize )
				return ( ArrayP::operator()(i,j)>>offJ );
			else
				return ( (ArrayP::operator()(i,j) & (p2[realN+offJ]-1))>>offJ );
		else // !offJ && realN < packetSize
			return ( ArrayP::operator()(i,j) & (p2[realN]-1) );
	else // j==ArrayP::n()-1 && moduloB(realN+offJ)
		return ( ArrayP::operator()(i,j) & (p2[moduloB(realN+offJ)]-1) );
	return 0;
}


/*inline*/ char binArray::block (int j) const
{	return binArray::block(0,j); 
}

/*inline*/ binArray& binArray::resize (const binArray& a)
{	
  ArrayP::resize(a); 
  realN=a.realN; 
  offJ=a.offJ; 

  return *this; 
}


//////////////////////
// Arrays modifiers //
//////////////////////

// simple assignment

binArray& binArray::operator  = (const boolean b) 
{	if ( !n() || !m() ) return *this;
	if ( !x() )	throw nullArray("binArray& binArray::operator=(const boolean)");
	packet* s; packet* t;
	packet* bom = x();
	packet* eol = x() + ArrayP::n();
	packet* eom = x() + (m()*l());
	register packet mask;
	if ( offJ && (offJ+n() < packetSize) )
	{	if ( b )
		{	mask = p2[offJ+n()] - p2[offJ];
			for (s=bom; s<eom; s+=l()) *s |= mask;
		}
		else
		{	mask = p2c[offJ+n()] + p2[offJ];
			for (s=bom; s<eom; s+=l()) *s &= mask;
		}
		return *this;
	}
	if ( offJ )
	{	if ( b )
		{	mask = p2c[offJ] + 1;
			for (s=bom; s<eom; s+=l()) *s |= mask;
		}
		else
		{	mask = p2[offJ] - 1;
			for (s=bom; s<eom; s+=l()) *s &= mask;
		}
		bom++;
	}
	if ( moduloB(offJ+n()) )
	{	eol--;
		if ( b )
		{	mask = p2[moduloB(offJ+n())] - 1;
			for (s=eol; s<eom; s+=l()) *s |= mask;
		}
		else
		{	mask = p2c[moduloB(offJ+n())] + 1;
			for (s=eol; s<eom; s+=l()) *s &= mask;
		}
	}
	mask = b ? fullP : 0;
	for (s=bom; s<eom; s+=l(), eol+=l())
		for (t=s; t<eol; t++) *t = mask;
	return *this;
}



binArray& binArray::operator= (const binArray& a)
{	
  if ( n()!=a.n() || m()!=a.m() ) 
    throw whereDimMismatch("binArray& binArray::operator=(const binArray&)",
			   m(),n(),a.m(),a.n());
  if ( !n() || !m() ) return *this;
  if ( !a.x() || !x() ) 
    throw nullArray("binArray& binArray::operator=(const binArray&)");

  if ( offJ != a.offJ )
    {	
      for (int i=0; i<m(); i++)
	for (int j=0; j<n(); j++)
	  set(i,j,a(i,j));
      return *this;
    }

  // n()==a.n() && offJ==a.offJ
  packet* s; packet* t; packet* ms; packet *mt;
  packet* bom = x();					
  packet* mbom = a.x();
  packet* eol = x() + ArrayP::n();	
  packet* meol = a.x() + a.ArrayP::n();
  packet* eom = x() + (m()*l());
  if ( offJ && (offJ+n() < packetSize) )
    {	
      register packet mask = p2[offJ+n()] - p2[offJ];
      for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l()) 
	{ 
	  *s &= ~mask; 
	  *s |= (*ms & mask); 
	}
      return *this;
    }

  if ( offJ )
    {	
      register packet mask = p2c[offJ] + 1; 
      for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l()) 
	{ 
	  *s &= ~mask; 
	  *s |= (*ms & mask); 
	}
      bom++; mbom++;
    }
  if ( moduloB(offJ+n()) )
    {	
      eol--; meol--;
      register packet mask = p2[moduloB(offJ+n())] - 1;
      for (s=eol, ms=meol; s<eom; s+=l(), ms+=a.l()) 
	{ 
	  *s &= ~mask; 
	  *s |= (*ms & mask); 
	}
    }
  for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l(), eol+=l())
    for (t=s, mt=ms; t<eol; t++, mt++) 
      *t = *mt;


  return *this;
}



binArray& binArray::operator &= (const binArray& a)
{	if ( n()!=a.n() || m()!=a.m() ) 
	 	throw whereDimMismatch(
			"binArray& binArray::operator&=(const binArray&)",
			m(),n(),a.m(),a.n());
	if ( !n() || !m() ) return *this;
	if ( !a.x() || !x() ) 
		throw nullArray("binArray& binArray::operator&=(const binArray&)");
	if ( offJ != a.offJ )
	{	for (int i=0; i<m(); i++)
			for (int j=0; j<n(); j++)
				set( i, j, operator()(i,j) & a(i,j) );
		return *this;
	}
	// n()==a.n() && offJ==a.offJ
	packet* s; packet* t; packet* ms; packet *mt;
	packet* bom = x();					packet* mbom = a.x();
	packet* eol = x() + ArrayP::n();	packet* meol = a.x() + a.ArrayP::n();
	packet* eom = x() + (m()*l());
	if ( offJ && (offJ+n() < packetSize) )
	{	register packet mask = p2[offJ+n()] - p2[offJ];
		for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l()) 
			*s &= (*ms | mask);
		return *this;
	}
	if ( offJ )
	{	register packet mask = p2[offJ] - 1; 
		for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l()) 
			*s &= (*ms | mask);
		bom++; mbom++;
	}
	if ( moduloB(offJ+n()) )
	{	eol--; meol--;
		register packet mask = p2c[moduloB(offJ+n())] + 1;
		for (s=eol, ms=meol; s<eom; s+=l(), ms+=a.l()) 
			*s &= (*ms | mask);
	}
	for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l(), eol+=l())
		for (t=s, mt=ms; t<eol; t++, mt++) *t &= *mt;
	return *this;
}


binArray& binArray::operator |= (const binArray& a)
{	if ( n()!=a.n() || m()!=a.m() ) 
	 	throw whereDimMismatch(
			"binArray& binArray::operator|=(const binArray&)",
			m(),n(),a.m(),a.n());
	if ( !n() || !m() ) return *this;
	if ( !a.x() || !x() ) 
		throw nullArray("binArray& binArray::operator|=(const binArray&)");
	if ( offJ != a.offJ )
	{	for (int i=0; i<m(); i++)
			for (int j=0; j<n(); j++)
				set( i, j, operator()(i,j) | a(i,j));
		return *this;
	}
	// n()==a.n() && offJ==a.offJ
	packet* s; packet* t; packet* ms; packet *mt;
	packet* bom = x();					packet* mbom = a.x();
	packet* eol = x() + ArrayP::n();	packet* meol = a.x() + a.ArrayP::n();
	packet* eom = x() + (m()*l());
	if ( offJ && (offJ+n() < packetSize) )
	{	register packet mask = p2[offJ+n()] - p2[offJ];
		for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l()) 
			*s |= (*ms & mask);
		return *this;
	}
	if ( offJ )
	{	register packet mask = p2c[offJ] + 1; 
		for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l()) 
			*s |= (*ms & mask); 
		bom++; mbom++;
	}
	if ( moduloB(offJ+n()) )
	{	eol--; meol--;
		register packet mask = p2[moduloB(offJ+n())] - 1;
		for (s=eol, ms=meol; s<eom; s+=l(), ms+=a.l()) 
			*s |= (*ms & mask);
	}
	for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l(), eol+=l())
		for (t=s, mt=ms; t<eol; t++, mt++) *t |= *mt;
	return *this;
}


binArray& binArray::operator ^= (const binArray& a)
{	if ( n()!=a.n() || m()!=a.m() ) 
	 	throw whereDimMismatch(
			"binArray& binArray::operator^=(const binArray&)",
			m(),n(),a.m(),a.n());
	if ( !n() || !m() ) return *this;
	if ( !a.x() || !x() ) 
		throw nullArray("binArray& binArray::operator^=(const binArray&)");
	if ( offJ != a.offJ )
	{	for (int i=0; i<m(); i++)
			for (int j=0; j<n(); j++)
				set( i, j, operator()(i,j) ^ a(i,j));
		return *this;
	}
	// n()==a.n() && offJ==a.offJ
	packet* s; packet* t; packet* ms; packet *mt;
	packet* bom = x();					packet* mbom = a.x();
	packet* eol = x() + ArrayP::n();	packet* meol = a.x() + a.ArrayP::n();
	packet* eom = x() + (m()*l());
	if ( offJ && (offJ+n() < packetSize) )
	{	register packet mask = p2[offJ+n()] - p2[offJ];
		for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l()) 
			*s ^= (*ms & mask);
		return *this;
	}
	if ( offJ )
	{	register packet mask = p2c[offJ] + 1; 
		for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l()) 
			*s ^= (*ms & mask); 
		bom++; mbom++;
	}
	if ( moduloB(offJ+n()) )
	{	eol--; meol--;
		register packet mask = p2[moduloB(offJ+n())] - 1;
		for (s=eol, ms=meol; s<eom; s+=l(), ms+=a.l()) 
			*s ^= (*ms & mask);
	}
	for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l(), eol+=l())
		for (t=s, mt=ms; t<eol; t++, mt++) *t ^= *mt;
	return *this;
}


/*inline*/ binArray& binArray::setOn (int i, int j)
{	ArrayP::operator()(i,floorDB(j+offJ)) |=  p2[moduloB(j+offJ)]; 
	return *this;
}


/*inline*/ binArray& binArray::setOff (int i, int j)
{	ArrayP::operator()(i,floorDB(j+offJ)) &= ~p2[moduloB(j+offJ)]; 
	return *this;
}


/*inline*/ binArray& binArray::flip (int i, int j)
{	ArrayP::operator()(i,floorDB(j+offJ)) ^=  p2[moduloB(j+offJ)]; 
	return *this;
}


/**inline**/ binArray& binArray::set (int i, int j, boolean a)
{	return (a ? setOn(i,j) : setOff(i,j)); 
}


/*inline*/ binArray& binArray::setOn (int j)
{	ArrayP::operator()(floorDB(j+offJ)) |=  p2[moduloB(j+offJ)]; 
	return *this;
}


/*inline*/ binArray& binArray::setOff (int j)
{	ArrayP::operator()(floorDB(j+offJ)) &= ~p2[moduloB(j+offJ)]; 
	return *this;
}


/*inline*/ binArray& binArray::flip (int j)
{	ArrayP::operator()(floorDB(j+offJ)) ^=  p2[moduloB(j+offJ)]; 
	return *this;
}


/*inline*/ binArray& binArray::set (int j, boolean a)
{	return (a) ? setOn(j) : setOff(j);
}



////////////////////////
// Arrays comparators //
////////////////////////


boolean binArray::operator  < (const binArray& a) const
{	if ( n()!=a.n() || m()!=a.m() ) 
	 	throw whereDimMismatch(
			"boolean binArray::operator<(const binArray&) const",
			m(),n(),a.m(),a.n());
	if ( !n() || !m() ) return false;
	if ( !a.x() || !x() ) 
		throw nullArray("boolean binArray::operator<(const binArray&) const");
	boolean equal = true;
	if ( offJ != a.offJ )
	{	for (int i=0; i<m(); i++)
			for (int j=0; j<n(); j++)
				switch ( tcomp( operator()(i,j), a(i,j) ) ) 
				{	case -1: return false;
					case +1: equal = false;
				}
		return !equal;
	}
	// n()==a.n() && offJ==a.offJ
	packet* s; packet* t; packet* ms; packet *mt;
	packet* bom = x();					packet* mbom = a.x();
	packet* eol = x() + ArrayP::n();	packet* meol = a.x() + a.ArrayP::n();
	packet* eom = x() + (m()*l());
	if ( offJ )
	{	register packet mask = p2c[offJ] + 1; 
		for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l()) 
			if ( *s & ~*ms & mask ) return false;
			else if ( equal ) equal = !( ~*s & *ms & mask );
		bom++; mbom++;
	}
	if ( moduloB(offJ+n()) )
	{	eol--; meol--;
		register packet mask = p2[moduloB(offJ+n())] - 1;
		for (s=eol, ms=meol; s<eom; s+=l(), ms+=a.l()) 
			if ( *s & ~*ms & mask ) return false;
			else if ( equal ) equal = !( ~*s & *ms & mask );
	}
	for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l(), eol+=l())
		for (t=s, mt=ms; t<eol; t++, mt++) 
			if ( *t & ~*mt ) return false;
			else if ( equal ) equal = ( *t == *mt );
	return !equal;
}


boolean binArray::operator <= (const binArray& a) const
{	if ( n()!=a.n() || m()!=a.m() ) 
	 	throw whereDimMismatch(
			"boolean binArray::operator<=(const binArray&) const",
			m(),n(),a.m(),a.n());
	if ( !n() || !m() ) return true;
	if ( !a.x() || !x() ) 
		throw nullArray("boolean binArray::operator<=(const binArray&) const");
	if ( offJ != a.offJ )
	{	for (int i=0; i<m(); i++)
			for (int j=0; j<n(); j++)
				if ( operator()(i,j) > a(i,j) ) return false;
		return true;
	}
	// n()==a.n() && offJ==a.offJ
	packet* s; packet* t; packet* ms; packet *mt;
	packet* bom = x();					packet* mbom = a.x();
	packet* eol = x() + ArrayP::n();	packet* meol = a.x() + a.ArrayP::n();
	packet* eom = x() + (m()*l());
	if ( offJ )
	{	register packet mask = p2c[offJ] + 1; 
		for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l()) 
			if ( *s & ~*ms & mask ) return false;
		bom++; mbom++;
	}
	if ( moduloB(offJ+n()) )
	{	eol--; meol--;
		register packet mask = p2[moduloB(offJ+n())] - 1;
		for (s=eol, ms=meol; s<eom; s+=l(), ms+=a.l()) 
			if ( *s & ~*ms & mask ) return false;
	}
	for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l(), eol+=l())
		for (t=s, mt=ms; t<eol; t++, mt++) 
			if ( *t & ~*mt ) return false;
	return true;
}


/*inline*/ boolean binArray::operator >  (const binArray& a) const
{	return (a < *this); 
}


/*inline*/ boolean binArray::operator >= (const binArray& a) const
{	return (a <= *this); 
}


boolean binArray::operator == (const binArray& a) const
{	if ( n()!=a.n() || m()!=a.m() ) 
	 	throw whereDimMismatch(
			"boolean binArray::operator==(const binArray&) const",
			m(),n(),a.m(),a.n());
	if ( !n() || !m() ) return true;
	if ( !a.x() || !x() ) 
		throw nullArray("boolean binArray::operator==(const binArray&) const");
	if ( offJ != a.offJ )
	{	for (int i=0; i<m(); i++)
			for (int j=0; j<n(); j++)
				if ( operator()(i,j) != a(i,j) ) return false;
		return true;
	}
	// n()==a.n() && offJ==a.offJ
	packet* s; packet* t; packet* ms; packet *mt;
	packet* bom = x();					packet* mbom = a.x();
	packet* eol = x() + ArrayP::n();	packet* meol = a.x() + a.ArrayP::n();
	packet* eom = x() + (m()*l());
	if ( offJ )
	{	register packet mask = p2c[offJ] + 1; 
		for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l()) 
			if ( (*s & mask) != (*ms & mask) ) return false;
		bom++; mbom++;
	}
	if ( moduloB(offJ+n()) )
	{	eol--; meol--;
		register packet mask = p2[moduloB(offJ+n())] - 1;
		for (s=eol, ms=meol; s<eom; s+=l(), ms+=a.l()) 
			if ( (*s & mask) != (*ms & mask) ) return false;
	}
	for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l(), eol+=l())
		for (t=s, mt=ms; t<eol; t++, mt++) 
			if ( *t != *mt ) return false;
	return true;
}


/*inline*/ boolean binArray::operator != (const binArray& a) const
{	return !(*this == a); 
}


int binArray::operator %= (const binArray& a) const
{	if ( n()!=a.n() || m()!=a.m() ) 
	 	throw whereDimMismatch(
			"int binArray::operator%=(const binArray&) const ",
			m(),n(),a.m(),a.n());
	if ( !n() || !m() ) return 0;
	if ( !a.x() || !x() ) 
		throw nullArray("int binArray::operator%=(const binArray&) const");
	int result = 0;
	int curcomp;
	if ( offJ != a.offJ )
	{	for (int i=0; i<m(); i++)
			for (int j=0; j<n(); j++)
			{	curcomp = tcomp( operator()(i,j), a(i,j) );
				if ( curcomp != result )
					if ( !result ) result = curcomp;
					else if ( curcomp ) return 2;
			}
		return result;
	}
	// n()==a.n() && offJ==a.offJ
	packet* s; packet* t; packet* ms; packet *mt;
	packet* bom = x();					packet* mbom = a.x();
	packet* eol = x() + ArrayP::n();	packet* meol = a.x() + a.ArrayP::n();
	packet* eom = x() + (m()*l());
	if ( offJ )
	{	register packet mask = p2c[offJ] + 1; 
		for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l()) 
		{	packet si = *s & mask;		packet msi = *ms & mask;
			if ( si != msi )
				if ( !(si & ~msi) ) // si < msi
					if ( result < 0 ) return 2;
					else result = +1;
				else if ( !(~si & msi) ) // si > msi
					if ( result > 0 ) return 2;
					else result = -1;
				else return 2;
		}
		bom++; mbom++;
	}
	if ( moduloB(offJ+n()) )
	{	eol--; meol--;
		register packet mask = p2[moduloB(offJ+n())] - 1;
		for (s=eol, ms=meol; s<eom; s+=l(), ms+=a.l())
		{	packet si = *s & mask;		packet msi = *ms & mask;
			if ( si != msi )
				if ( !(si & ~msi) ) // si < msi
					if ( result < 0 ) return 2;
					else result = +1;
				else if ( !(~si & msi) ) // si > msi
					if ( result > 0 ) return 2;
					else result = -1;
				else return 2;
		}
	}
	for (s=bom, ms=mbom; s<eom; s+=l(), ms+=a.l(), eol+=l())
		for (t=s, mt=ms; t<eol; t++, mt++) 
			if ( *t != *mt )
				if ( !(*t & ~*mt) ) // *t < *mt
					if ( result < 0 ) return 2;
					else result = +1;
				else if ( !(~*t & *mt) ) // *t > *mt
					if ( result > 0 ) return 2;
					else result = -1;
				else return 2;
	return result;
}



//////////////////////
// Arrays operators //
//////////////////////

binMatrix binArray::operator & (const binArray& a) const
{	if ( n()!=a.n() || m()!=a.m() ) 
	 	throw whereDimMismatch(
			"binMatrix binArray::operator&(const binArray&) const",
			m(),n(),a.m(),a.n());
	if ( !n() || !m() ) return binMatrix(m(),n());
	if ( !a.x() || !x() ) 
		throw nullArray("binMatrix binArray::operator&(const binArray&) const");
	if ( offJ != a.offJ )
	{	binMatrix result(m(),n());
		for (int i=0; i<m(); i++)
			for (int j=0; j<n(); j++)
				result.set( i, j, operator()(i,j) && a(i,j) );
		return result;
	}
	// n()==a.n() && offJ==a.offJ
	binMatrix result(binArray(m(),n(),offJ));
  	packet* s; packet* t; packet* ms; packet* mt; packet* rt;
	packet* eol = x() + ArrayP::n();	
	packet* eom = x() + (m()*l());
	for (s=x(), ms=a.x(), rt=result.x(); s<eom; s+=l(), ms+=a.l(), eol+=l())
		for (t=s, mt=ms; t<eol; t++, rt++, mt++) 
			*rt = (*t & *mt);
	return result;
}


binMatrix binArray::operator | (const binArray& a) const
{	if ( n()!=a.n() || m()!=a.m() ) 
	 	throw whereDimMismatch(
			"binMatrix binArray::operator|(const binArray&) const",
			m(),n(),a.m(),a.n());
	if ( !n() || !m() ) return binMatrix(m(),n());
	if ( !a.x() || !x() ) 
		throw nullArray("binMatrix binArray::operator&(const binArray&) const");
	if ( offJ != a.offJ )
	{	binMatrix result(m(),n());
		for (int i=0; i<m(); i++)
			for (int j=0; j<n(); j++)
				result.set( i, j, operator()(i,j) || a(i,j) );
		return result;
	}
	// n()==a.n() && offJ==a.offJ
	binMatrix result(binArray(m(),n(),offJ));
	packet* s; packet* t; packet* ms; packet* mt; packet* rt;
	packet* eol = x() + ArrayP::n();	
	packet* eom = x() + (m()*l());
	for (s=x(), ms=a.x(), rt=result.x(); s<eom; s+=l(), ms+=a.l(), eol+=l())
		for (t=s, mt=ms; t<eol; t++, rt++, mt++) 
			*rt = (*t | *mt);
	return result;
}


binMatrix binArray::operator ^ (const binArray& a) const
{	if ( n()!=a.n() || m()!=a.m() ) 
	 	throw whereDimMismatch(
			"binMatrix binArray::operator^(const binArray&) const",
			m(),n(),a.m(),a.n());
	if ( !n() || !m() ) return binMatrix(m(),n());
	if ( !a.x() || !x() ) 
		throw nullArray("binMatrix binArray::operator^(const binArray&) const");
	if ( offJ != a.offJ )
	{	binMatrix result(m(),n());
		for (int i=0; i<m(); i++)
			for (int j=0; j<n(); j++)
				result.set( i, j, ( operator()(i,j) != a(i,j) ) );
		return result;
	}
	// n()==a.n() && offJ==a.offJ
	binMatrix result(binArray(m(),n(),offJ));
	packet* s; packet* t; packet* ms; packet* mt; packet* rt;
	packet* eol = x() + ArrayP::n();	
	packet* eom = x() + (m()*l());
	for (s=x(), ms=a.x(), rt=result.x(); s<eom; s+=l(), ms+=a.l(), eol+=l())
		for (t=s, mt=ms; t<eol; t++, rt++, mt++) 
			*rt = (*t ^ *mt);
	return result;
}


binMatrix binArray::operator ~ () const
{	if ( !n() || !m() ) return binMatrix(m(),n());
	if ( !x() ) throw nullArray("binMatrix binArray::operator~() const");
	binMatrix result(binArray(m(),n(),offJ));
	packet* s; packet* t; packet* rt;
	packet* eol = x() + ArrayP::n();	
	packet* eom = x() + (m()*l());
	for (s=x(), rt=result.x(); s<eom; s+=l(), eol+=l())
		for (t=s; t<eol; t++, rt++) 
			*rt = ~(*t);
	return result;
}



binMatrix binArray::andByRows () const
{	if ( !n() ) return binMatrix(1,0);
	if ( !m() ) return binMatrix(1,n(),true);
	if ( !x() ) throw nullArray("binMatrix binArray::andByRows() const");
	binMatrix result(binArray(1,n(),offJ));
	packet* s; packet* t; packet* rt;
	packet* eol = x() + ArrayP::n();
	packet* eom = x() + (m()*l());
	for (t=x(), rt=result.x(); t<eol; t++, rt++)
	{	*rt = fullP;
		for (s=t; s<eom; s+=l())
			if ( !(*rt = (*rt & *s)) ) break;
	}
	return result;
}


binMatrix binArray::orByRows () const
{	if ( !n() ) return binMatrix(1,0);
	if ( !m() ) return binMatrix(1,n(),false);
	if ( !x() ) throw nullArray("binMatrix binArray::orByRows() const");
	binMatrix result(binArray(1,n(),offJ));
	packet* s; packet* t; packet* rt;
	packet* eol = x() + ArrayP::n();
	packet* eom = x() + (m()*l());
	for (t=x(), rt=result.x(); t<eol; t++, rt++)
	{	*rt = 0;
		for (s=t; s<eom; s+=l())
			if ( (*rt = (*rt | *s)) == fullP ) break;
	}
	return result;
}


binMatrix binArray::xorByRows () const
{	if ( !n() ) return binMatrix(1,0);
	if ( !m() ) return binMatrix(1,n(),false);
	if ( !x() ) throw nullArray("binMatrix binArray::xorByRows() const");
	binMatrix result(binArray(1,n(),offJ));
	packet* s; packet* t; packet* rt;
	packet* eol = x() + ArrayP::n();
	packet* eom = x() + (m()*l());
	for (t=x(), rt=result.x(); t<eol; t++, rt++)
	{	*rt = 0;
		for (s=t; s<eom; s+=l())
			*rt = *rt ^ *s;
	}
	return result;
}


binMatrix binArray::andByCols () const
{	if ( !m() ) return binMatrix(1,0);
	if ( !n() ) return binMatrix(1,m(),true);
	if ( !x() ) throw nullArray("binMatrix binArray::andByCols() const");
	binMatrix result(1,m(),true);
	packet* s; packet* t; int rt;
	packet* eol = x() + ArrayP::n() - 1;
	packet* eom = x() + (m()*l());
	packet maskL= p2c[offJ] + 1;
	packet maskR= (moduloB(offJ+n()) ? p2[moduloB(offJ+n())] - 1 : fullP);
	for (s=x(), rt=0; s<eom; s+=l(), rt++, eol+=l())
		if ( *s < maskL ) result.setOff(rt);
		else if ( eol > s ) 
			if  ( (*eol & maskR) < maskR ) result.setOff(rt);
			else 
				for (t=s+1; t<eol; t++)
					if ( *t < fullP ) { result.setOff(rt); break; }
	return result;
}


binMatrix binArray::orByCols () const
{	if ( !m() ) return binMatrix(1,0);
	if ( !n() ) return binMatrix(1,m(),false);
	if ( !x() ) throw nullArray("binMatrix binArray::orByCols() const");
	binMatrix result(1,m(),false);
	packet* s; packet* t; int rt;
	packet* eol = x() + ArrayP::n() - 1;
	packet* eom = x() + (m()*l());
	packet maskL= p2c[offJ] + 1;
	packet maskR= (moduloB(offJ+n()) ? p2[moduloB(offJ+n())] - 1 : fullP);
	for (s=x(), rt=0; s<eom; s+=l(), rt++, eol+=l())
		if ( *s & maskL ) result.setOn(rt);
		else if ( eol > s )
			if ( *eol & maskR ) result.setOn(rt);
			else 
				for (t=s+1; t<eol; t++)
					if ( *t ) { result.setOn(rt); break; }
	return result;
}


binMatrix binArray::xorByCols () const
{	if ( !m() ) return binMatrix(1,0);
	binMatrix result(1,m(),false);
	for (int i=0; i<m(); i++)
		for (int j=0; j<n(); j++)
			if ( operator()(i,j) ) result.flip(i);
	return result;
}


binMatrix binArray::t () const
{	binMatrix result(n(),m());
	for (int i=0; i<m(); i++)
		for (int j=0; j<n(); j++)
			result.set(j,i,operator()(i,j));
	return result;
}


///////////////////
// I/O of arrays //
///////////////////

ostream& operator << (ostream& s, const binArray& a)
{	for (int i=0; i<a.m(); i++)
	{	for (int j=0; j<a.n(); j++)
			s << ( a(i,j) ? "1" : "0" );
		s << endl;
	}
	return s;
}

istream& operator >> (istream& s, binArray& a)
{	char t;
	for (int i=0; i<a.m(); i++)
		for (int j=0; j<a.n(); j++)
		{	while ( s.get(t) && t!='0' && t!='1' ) ;
			if ( ( !s && (i<a.m()-1 || j<a.n()-1) ) || ( t!='0' && t!='1' ) )
			{	cerr<< "\n### Bad format in istream& operator>>(istream&,"
					<< "binArray&) !\n";
				return s;
			}
			if ( t=='1' ) a.setOn(i,j); else a.setOff(i,j);
		}
	return s;
}


///////////////////////////
// Matrices constructors //
///////////////////////////

/*inline*/ binMatrix::binMatrix (int n) :
	ArrayP(ceilDB(n), 1, ceilDB(n), new packet [ceilDB(n)]), binArray(1, n) { 
}

/**inline**/ binMatrix::binMatrix (int m,	int n) :
	ArrayP(ceilDB(n), m, ceilDB(n), new packet [m * ceilDB(n)]),
	binArray(m, n) { 
}

binMatrix::binMatrix (int m, int n, const boolean b) :
	ArrayP(ceilDB(n), m, ceilDB(n), new packet [m * ceilDB(n)]),
	binArray(m, n) 
{	
  binArray::operator=(b);
}

binMatrix::binMatrix (int n, const boolean b) :
	ArrayP(ceilDB(n), 1, ceilDB(n), new packet [ceilDB(n)]),
	binArray(1, n) 
{	binArray::operator=(b); 
}

binMatrix::binMatrix (const binMatrix& a) :
	ArrayP(ceilDB(a.n()),a.m(),ceilDB(a.n()), new packet[a.m()*ceilDB(a.n())]),
	binArray(a.m(), a.n(), 0)
{	binArray::operator=(a); 
}

binMatrix::binMatrix (const binArray& a) :
	ArrayP(ceilDB(a.n()),a.m(),ceilDB(a.n()), new packet[a.m()*ceilDB(a.n())]),
	binArray(a.m(), a.n(), 0)
{	try { binArray::operator=(a); }
	catch (nullArray) {}
}

/**inline**/ binMatrix::binMatrix () : binArray() { 
}

/*inline*/ binMatrix::~binMatrix ()
{	delete [] x(); 
}



////////////////////////
// Matrices modifiers //
////////////////////////

binMatrix& binMatrix::resize (const binArray& a)
{	
  delete [] this->x(); 
  binArray::resize(a); 
  this->x(new packet [a.m()*a.n()]);

  //*this = a;
  operator= (a);  


  return *this; 
}


binMatrix& binMatrix::resize (int m, int n, const boolean b)
{	MatrixP::resize(m, ceilDB(n));
 	packet  mask = b ? fullP : 0;
	packet* eom = x()+(m*ceilDB(n));
	for (packet* t=x(); t<eom; t++) *t = mask;
	realN = n; offJ=0;
	return *this;
}


binMatrix& binMatrix::resize (int n, const boolean b)
{	MatrixP::resize(1, ceilDB(n));
 	packet  mask = b ? fullP : 0;
	packet* eom = x()+ceilDB(n);
	for (packet* t=x(); t<eom; t++) *t = mask;
	realN = n; offJ=0;
	return *this;
}


binMatrix& binMatrix::resize (int m, int n)
{	MatrixP::resize(m, ceilDB(n));
	realN = n; offJ=0;
	return *this;
}


binMatrix& binMatrix::resize (int n/**=0**/)
{	MatrixP::resize(1, ceilDB(n));
	realN = n; offJ=0;
	return *this;
}


binMatrix& binMatrix::resize_(int m, int n)
{	MatrixP::resize_(m, ceilDB(n+offJ));
	realN = n;
	return *this;
}

binMatrix& binMatrix::resize_(int n)
{	MatrixP::resize_(ceilDB(n+offJ));
	realN = n;
	return *this;
}



////////////////////////
// Matrices operators //
////////////////////////

binMatrix& binMatrix::operator = (const binMatrix& a)
{	binArray::operator=(a); return *this; 
}

binMatrix& binMatrix::operator = (const binArray& a)
{	binArray::operator=(a); return *this; 
}

binMatrix& binMatrix::operator = (const boolean a)
{	binArray::operator=(a); return *this; 
}