for (i=0; i<smallerSize; i++)
		if (reg[i] != rhs.reg[i]) return false;

	for (i=smallerSize; i<reg.size(); i++)
		if (reg[i] != 0) return false;

	for (i=smallerSize; i<rhs.reg.size(); i++)
		if (rhs.reg[i] != 0) return false;

	return true;
}

std::ostream& operator<<(std::ostream& out, const PolynomialMod2 &a)
{
	// Get relevant conversion specifications from ostream.
	long f = out.flags() & std::ios::basefield;	// Get base digits.
	int bits, block;
	char suffix;
	switch(f)
	{
	case std::ios::oct :
		bits = 3;
		block = 4;
		suffix = 'o';
		break;
	case std::ios::hex :
		bits = 4;
		block = 2;
		suffix = 'h';
		break;
	default :
		bits = 1;
		block = 8;
		suffix = 'b';
	}

	if (!a)
		return out << '0' << suffix;

	SecBlock<char> s(a.BitCount()/bits+1);
	unsigned i;
	const char vec[]="0123456789ABCDEF";

	for (i=0; i*bits < a.BitCount(); i++)
	{
		int digit=0;
		for (int j=0; j<bits; j++)
			digit |= a[i*bits+j] << j;
		s[i]=vec[digit];
	}

	while (i--)
	{
		out << s[i];
		if (i && (i%block)==0)
			out << ',';
	}

	return out << suffix;
}

PolynomialMod2 PolynomialMod2::Gcd(const PolynomialMod2 &a, const PolynomialMod2 &b)
{
	return EuclideanDomainOf<PolynomialMod2>().Gcd(a, b);
}

PolynomialMod2 PolynomialMod2::InverseMod(const PolynomialMod2 &modulus) const
{
	typedef EuclideanDomainOf<PolynomialMod2> Domain;
	return QuotientRing<Domain>(Domain(), modulus).MultiplicativeInverse(*this);
}

bool PolynomialMod2::IsIrreducible() const
{
	signed int d = Degree();
	if (d <= 0)
		return false;

	PolynomialMod2 t(2), u(t);
	for (int i=1; i<=d/2; i++)
	{
		u = u.Squared()%(*this);
		if (!Gcd(u+t, *this).IsUnit())
			return false;
	}
	return true;
}

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

GF2NP::GF2NP(const PolynomialMod2 &modulus)
	: QuotientRing<EuclideanDomainOf<PolynomialMod2> >(EuclideanDomainOf<PolynomialMod2>(), modulus), m(modulus.Degree()) 
{
}

GF2NP::Element GF2NP::SquareRoot(const Element &a) const
{
	Element r = a;
	for (unsigned int i=1; i<m; i++)
		r = Square(r);
	return r;
}

GF2NP::Element GF2NP::HalfTrace(const Element &a) const
{
	assert(m%2 == 1);
	Element h = a;
	for (unsigned int i=1; i<=(m-1)/2; i++)
		h = Add(Square(Square(h)), a);
	return h;
}

GF2NP::Element GF2NP::SolveQuadraticEquation(const Element &a) const
{
	if (m%2 == 0)
	{
		Element z, w;
		do
		{
			LC_RNG rng(11111);
			Element p(rng, m);
			z = PolynomialMod2::Zero();
			w = p;
			for (unsigned int i=1; i<=m-1; i++)
			{
				w = Square(w);
				z = Square(z);
				Accumulate(z, Multiply(w, a));
				Accumulate(w, p);
			}
		} while (w.IsZero());
		return z;
	}
	else
		return HalfTrace(a);
}

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

GF2NT::GF2NT(unsigned int t0, unsigned int t1, unsigned int t2)
	: GF2NP(PolynomialMod2::Trinomial(t0, t1, t2))
	, t0(t0), t1(t1)
	, result((word)0, m)
{
	assert(t0 > t1 && t1 > t2 && t2==0);
}

const GF2NT::Element& GF2NT::MultiplicativeInverse(const Element &a) const
{
	if (t0-t1 < WORD_BITS)
		return GF2NP::MultiplicativeInverse(a);

	SecWordBlock T(m_modulus.reg.size() * 4);
	word *b = T;
	word *c = T+m_modulus.reg.size();
	word *f = T+2*m_modulus.reg.size();
	word *g = T+3*m_modulus.reg.size();
	unsigned int bcLen=1, fgLen=m_modulus.reg.size();
	unsigned int k=0;

	SetWords(T, 0, 3*m_modulus.reg.size());
	b[0]=1;
	assert(a.reg.size() <= m_modulus.reg.size());
	CopyWords(f, a.reg, a.reg.size());
	CopyWords(g, m_modulus.reg, m_modulus.reg.size());

	while (1)
	{
		word t=f[0];
		while (!t)
		{
			ShiftWordsRightByWords(f, fgLen, 1);
			if (c[bcLen-1])
				bcLen++;
			assert(bcLen <= m_modulus.reg.size());
			ShiftWordsLeftByWords(c, bcLen, 1);
			k+=WORD_BITS;
			t=f[0];
		}

		unsigned int i=0;
		while (t%2 == 0)
		{
			t>>=1;
			i++;
		}
		k+=i;

		if (t==1 && CountWords(f, fgLen)==1)
			break;

		if (i==1)
		{
			ShiftWordsRightByBits(f, fgLen, 1);
			t=ShiftWordsLeftByBits(c, bcLen, 1);
		}
		else
		{
			ShiftWordsRightByBits(f, fgLen, i);
			t=ShiftWordsLeftByBits(c, bcLen, i);
		}
		if (t)
		{
			c[bcLen] = t;
			bcLen++;
			assert(bcLen <= m_modulus.reg.size());
		}

		if (f[fgLen-1]==0 && g[fgLen-1]==0)
			fgLen--;

		if (f[fgLen-1] < g[fgLen-1])
		{
			std::swap(f, g);
			std::swap(b, c);
		}

		XorWords(f, g, fgLen);
		XorWords(b, c, bcLen);
	}

	while (k >= WORD_BITS)
	{
		word temp = b[0];
		// right shift b
		for (unsigned i=0; i+1<BitsToWords(m); i++)
			b[i] = b[i+1];
		b[BitsToWords(m)-1] = 0;

		if (t1 < WORD_BITS)
			for (unsigned int j=0; j<WORD_BITS-t1; j++)
				temp ^= ((temp >> j) & 1) << (t1 + j);
		else
			b[t1/WORD_BITS-1] ^= temp << t1%WORD_BITS;

		if (t1 % WORD_BITS)
			b[t1/WORD_BITS] ^= temp >> (WORD_BITS - t1%WORD_BITS);

		if (t0%WORD_BITS)
		{
			b[t0/WORD_BITS-1] ^= temp << t0%WORD_BITS;
			b[t0/WORD_BITS] ^= temp >> (WORD_BITS - t0%WORD_BITS);
		}
		else
			b[t0/WORD_BITS-1] ^= temp;

		k -= WORD_BITS;
	}

	if (k)
	{
		word temp = b[0] << (WORD_BITS - k);
		ShiftWordsRightByBits(b, BitsToWords(m), k);

		if (t1 < WORD_BITS)
			for (unsigned int j=0; j<WORD_BITS-t1; j++)
				temp ^= ((temp >> j) & 1) << (t1 + j);
		else
			b[t1/WORD_BITS-1] ^= temp << t1%WORD_BITS;

		if (t1 % WORD_BITS)
			b[t1/WORD_BITS] ^= temp >> (WORD_BITS - t1%WORD_BITS);

		if (t0%WORD_BITS)
		{
			b[t0/WORD_BITS-1] ^= temp << t0%WORD_BITS;
			b[t0/WORD_BITS] ^= temp >> (WORD_BITS - t0%WORD_BITS);
		}
		else
			b[t0/WORD_BITS-1] ^= temp;
	}

	CopyWords(result.reg.begin(), b, result.reg.size());
	return result;
}

const GF2NT::Element& GF2NT::Multiply(const Element &a, const Element &b) const
{
	unsigned int aSize = STDMIN(a.reg.size(), result.reg.size());
	Element r((word)0, m);

	for (int i=m-1; i>=0; i--)
	{
		if (r[m-1])
		{
			ShiftWordsLeftByBits(r.reg.begin(), r.reg.size(), 1);
			XorWords(r.reg.begin(), m_modulus.reg, r.reg.size());
		}
		else
			ShiftWordsLeftByBits(r.reg.begin(), r.reg.size(), 1);

		if (b[i])
			XorWords(r.reg.begin(), a.reg, aSize);
	}

	if (m%WORD_BITS)
		r.reg.begin()[r.reg.size()-1] = (word)Crop(r.reg[r.reg.size()-1], m%WORD_BITS);

	CopyWords(result.reg.begin(), r.reg.begin(), result.reg.size());
	return result;
}

const GF2NT::Element& GF2NT::Reduced(const Element &a) const
{
	if (t0-t1 < WORD_BITS)
		return m_domain.Mod(a, m_modulus);

	SecWordBlock b(a.reg);

	unsigned i;
	for (i=b.size()-1; i>=BitsToWords(t0); i--)
	{
		word temp = b[i];

		if (t0%WORD_BITS)
		{
			b[i-t0/WORD_BITS] ^= temp >> t0%WORD_BITS;
			b[i-t0/WORD_BITS-1] ^= temp << (WORD_BITS - t0%WORD_BITS);
		}
		else
			b[i-t0/WORD_BITS] ^= temp;

		if ((t0-t1)%WORD_BITS)
		{
			b[i-(t0-t1)/WORD_BITS] ^= temp >> (t0-t1)%WORD_BITS;
			b[i-(t0-t1)/WORD_BITS-1] ^= temp << (WORD_BITS - (t0-t1)%WORD_BITS);
		}
		else
			b[i-(t0-t1)/WORD_BITS] ^= temp;
	}

	if (i==BitsToWords(t0)-1 && t0%WORD_BITS)
	{
		word mask = ((word)1<<(t0%WORD_BITS))-1;
		word temp = b[i] & ~mask;
		b[i] &= mask;

		b[i-t0/WORD_BITS] ^= temp >> t0%WORD_BITS;

		if ((t0-t1)%WORD_BITS)
		{
			b[i-(t0-t1)/WORD_BITS] ^= temp >> (t0-t1)%WORD_BITS;
			if ((t0-t1)%WORD_BITS > t0%WORD_BITS)
				b[i-(t0-t1)/WORD_BITS-1] ^= temp << (WORD_BITS - (t0-t1)%WORD_BITS);
			else
				assert(temp << (WORD_BITS - (t0-t1)%WORD_BITS) == 0);
		}
		else
			b[i-(t0-t1)/WORD_BITS] ^= temp;
	}

	SetWords(result.reg.begin(), 0, result.reg.size());
	CopyWords(result.reg.begin(), b, STDMIN(b.size(), result.reg.size()));
	return result;
}

void GF2NP::DEREncodeElement(BufferedTransformation &out, const Element &a) const
{
	a.DEREncodeAsOctetString(out, MaxElementByteLength());
}

void GF2NP::BERDecodeElement(BufferedTransformation &in, Element &a) const
{
	a.BERDecodeAsOctetString(in, MaxElementByteLength());
}

void GF2NT::DEREncode(BufferedTransformation &bt) const
{
	DERSequenceEncoder seq(bt);
		ASN1::characteristic_two_field().DEREncode(seq);
		DERSequenceEncoder parameters(seq);
			DEREncodeUnsigned(parameters, m);
			ASN1::tpBasis().DEREncode(parameters);
			DEREncodeUnsigned(parameters, t1);
		parameters.MessageEnd();
	seq.MessageEnd();
}

void GF2NPP::DEREncode(BufferedTransformation &bt) const
{
	DERSequenceEncoder seq(bt);
		ASN1::characteristic_two_field().DEREncode(seq);
		DERSequenceEncoder parameters(seq);
			DEREncodeUnsigned(parameters, m);
			ASN1::ppBasis().DEREncode(parameters);
			DERSequenceEncoder pentanomial(parameters);
				DEREncodeUnsigned(pentanomial, t3);
				DEREncodeUnsigned(pentanomial, t2);
				DEREncodeUnsigned(pentanomial, t1);
			pentanomial.MessageEnd();
		parameters.MessageEnd();
	seq.MessageEnd();
}

GF2NP * BERDecodeGF2NP(BufferedTransformation &bt)
{
	// VC60 workaround: auto_ptr lacks reset()
	member_ptr<GF2NP> result;

	BERSequenceDecoder seq(bt);
		if (OID(seq) != ASN1::characteristic_two_field())
			BERDecodeError();
		BERSequenceDecoder parameters(seq);
			unsigned int m;
			BERDecodeUnsigned(parameters, m);
			OID oid(parameters);
			if (oid == ASN1::tpBasis())
			{
				unsigned int t1;
				BERDecodeUnsigned(parameters, t1);
				result.reset(new GF2NT(m, t1, 0));
			}
			else if (oid == ASN1::ppBasis())
			{
				unsigned int t1, t2, t3;
				BERSequenceDecoder pentanomial(parameters);
				BERDecodeUnsigned(pentanomial, t3);
				BERDecodeUnsigned(pentanomial, t2);
				BERDecodeUnsigned(pentanomial, t1);
				pentanomial.MessageEnd();
				result.reset(new GF2NPP(m, t3, t2, t1, 0));
			}
			else
			{
				BERDecodeError();
				return NULL;
			}
		parameters.MessageEnd();
	seq.MessageEnd();

	return result.release();
}

NAMESPACE_END