}
	return *this;
}

Integer& Integer::operator<<=(unsigned int n)
{
	const unsigned int wordCount = WordCount();
	const unsigned int shiftWords = n / WORD_BITS;
	const unsigned int shiftBits = n % WORD_BITS;

	reg.CleanGrow(RoundupSize(wordCount+bitsToWords(n)));
	ShiftWordsLeftByWords(reg, wordCount + shiftWords, shiftWords);
	ShiftWordsLeftByBits(reg+shiftWords, wordCount+bitsToWords(shiftBits), shiftBits);
	return *this;
}

Integer& Integer::operator>>=(unsigned int n)
{
	const unsigned int wordCount = WordCount();
	const unsigned int shiftWords = n / WORD_BITS;
	const unsigned int shiftBits = n % WORD_BITS;

	ShiftWordsRightByWords(reg, wordCount, shiftWords);
	if (wordCount > shiftWords)
		ShiftWordsRightByBits(reg, wordCount-shiftWords, shiftBits);
	if (IsNegative() && WordCount()==0)   // avoid -0
		*this = Zero();
	return *this;
}

void PositiveMultiply(Integer &product, const Integer &a, const Integer &b)
{
	unsigned aSize = RoundupSize(a.WordCount());
	unsigned bSize = RoundupSize(b.WordCount());

	product.reg.CleanNew(RoundupSize(aSize+bSize));
	product.sign = Integer::POSITIVE;

	SecWordBlock workspace(aSize + bSize);
	AsymmetricMultiply(product.reg, workspace, a.reg, aSize, b.reg, bSize);
}

void Multiply(Integer &product, const Integer &a, const Integer &b)
{
	PositiveMultiply(product, a, b);

	if (a.NotNegative() != b.NotNegative())
		product.Negate();
}

Integer operator*(const Integer &a, const Integer &b)
{
	Integer product;
	Multiply(product, a, b);
	return product;
}

/*
void PositiveDivide(Integer &remainder, Integer &quotient,
				   const Integer &dividend, const Integer &divisor)
{
	remainder.reg.CleanNew(divisor.reg.size);
	remainder.sign = Integer::POSITIVE;
	quotient.reg.New(0);
	quotient.sign = Integer::POSITIVE;
	unsigned i=dividend.BitCount();
	while (i--)
	{
		word overflow = ShiftWordsLeftByBits(remainder.reg, remainder.reg.size, 1);
		remainder.reg[0] |= dividend[i];
		if (overflow || remainder >= divisor)
		{
			Subtract(remainder.reg, remainder.reg, divisor.reg, remainder.reg.size);
			quotient.SetBit(i);
		}
	}
}
*/

void PositiveDivide(Integer &remainder, Integer &quotient,
				   const Integer &a, const Integer &b)
{
	unsigned aSize = a.WordCount();
	unsigned bSize = b.WordCount();

	if (!bSize)
		throw Integer::DivideByZero();

	if (a.PositiveCompare(b) == -1)
	{
		remainder = a;
		remainder.sign = Integer::POSITIVE;
		quotient = Integer::Zero();
		return;
	}

	aSize += aSize%2;   // round up to next even number
	bSize += bSize%2;

	remainder.reg.CleanNew(RoundupSize(bSize));
	remainder.sign = Integer::POSITIVE;
	quotient.reg.CleanNew(RoundupSize(aSize-bSize+2));
	quotient.sign = Integer::POSITIVE;

	SecWordBlock T(aSize+2*bSize+4);
	Divide(remainder.reg, quotient.reg, T, a.reg, aSize, b.reg, bSize);
}

void Integer::Divide(Integer &remainder, Integer &quotient, const Integer &dividend, const Integer &divisor)
{
	PositiveDivide(remainder, quotient, dividend, divisor);

	if (dividend.IsNegative())
	{
		quotient.Negate();
		if (!!remainder)
		{
			--quotient;
			remainder = divisor.AbsoluteValue() - remainder;
		}
	}

	if (divisor.IsNegative())
		quotient.Negate();
}

Integer operator/(const Integer &a, const Integer &b)
{
	Integer remainder, quotient;
	Integer::Divide(remainder, quotient, a, b);
	return quotient;
}

Integer operator%(const Integer &a, const Integer &b)
{
	Integer remainder, quotient;
	Integer::Divide(remainder, quotient, a, b);
	return remainder;
}

word Integer::ShortDivide(Integer &quotient, const Integer &dividend, word divisor)
{
	if (!divisor)
		throw Integer::DivideByZero();

	assert(divisor);

	if ((divisor & (divisor-1)) == 0)   // divisor is a power of 2
	{
		quotient = dividend >> (BitPrecision(divisor)-1);
		return dividend.reg[0] & (divisor-1);
	}

	unsigned int i = dividend.WordCount();
	quotient.reg.CleanNew(RoundupSize(i));
	word remainder = 0;
	while (i--)
	{
		quotient.reg[i] = word(MAKE_DWORD(dividend.reg[i], remainder) / divisor);
		remainder = word(MAKE_DWORD(dividend.reg[i], remainder) % divisor);
	}

	if (dividend.NotNegative())
		quotient.sign = POSITIVE;
	else
	{
		quotient.sign = NEGATIVE;
		if (remainder)
		{
			--quotient;
			remainder = divisor - remainder;
		}
	}

	return remainder;
}

Integer operator/(const Integer &a, word b)
{
	Integer quotient;
	Integer::ShortDivide(quotient, a, b);
	return quotient;
}

word operator%(const Integer &dividend, word divisor)
{
	if (!divisor)
		throw Integer::DivideByZero();

	assert(divisor);

	word remainder;

	if ((divisor & (divisor-1)) == 0)   // divisor is a power of 2
		remainder = dividend.reg[0] & (divisor-1);
	else
	{
		unsigned int i = dividend.WordCount();

		if (divisor <= 5)
		{
			dword sum=0;
			while (i--)
				sum += dividend.reg[i];
			remainder = word(sum%divisor);
		}
		else
		{
			remainder = 0;
			while (i--)
				remainder = word(MAKE_DWORD(dividend.reg[i], remainder) % divisor);
		}
	}

	if (dividend.IsNegative() && remainder)
		remainder = divisor - remainder;

	return remainder;
}

void Integer::Negate()
{
	if (!!(*this))  // don't flip sign if *this==0
		sign = Sign(1-sign);
}

int Integer::PositiveCompare(const Integer& t) const
{
	unsigned size = WordCount(), tSize = t.WordCount();

	if (size == tSize)
		return CryptoPP::Compare(reg, t.reg, size);
	else
		return size > tSize ? 1 : -1;
}

int Integer::Compare(const Integer& t) const
{
	if (NotNegative())
	{
		if (t.NotNegative())
			return PositiveCompare(t);
		else
			return 1;
	}
	else
	{
		if (t.NotNegative())
			return -1;
		else
			return -PositiveCompare(t);
	}
}

Integer Integer::SquareRoot() const
{
	if (!IsPositive())
		return Zero();

	// overestimate square root
	Integer x, y = Power2((BitCount()+1)/2);
	assert(y*y >= *this);

	do
	{
		x = y;
		y = (x + *this/x) >> 1;
	} while (y<x);

	return x;
}

bool Integer::IsSquare() const
{
	Integer r = SquareRoot();
	return *this == r.Squared();
}

bool Integer::IsUnit() const
{
	return (WordCount() == 1) && (reg[0] == 1);
}

Integer Integer::MultiplicativeInverse() const
{
	return IsUnit() ? *this : Zero();
}

Integer a_times_b_mod_c(const Integer &x, const Integer& y, const Integer& m)
{
	return x*y%m;
}

Integer a_exp_b_mod_c(const Integer &x, const Integer& e, const Integer& m)
{
	if (m.IsEven())
	{
		ModularArithmetic mr(m);
		return mr.ConvertOut(mr.Exponentiate(mr.ConvertIn(x), e));
	}
	else
	{
		MontgomeryRepresentation mr(m);
		return mr.ConvertOut(mr.Exponentiate(mr.ConvertIn(x), e));
	}
}

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

Integer Integer::InverseMod(const Integer &m) const
{
	assert(m.NotNegative());

	if (IsNegative() || *this>=m)
		return (*this%m).InverseMod(m);

	if (m.IsEven())
	{
		if (!m || IsEven())
			return Zero();	// no inverse
		if (*this == One())
			return One();

		Integer u = m.InverseMod(*this);
		return !u ? Zero() : (m*(*this-u)+1)/(*this);
	}

	SecBlock<word> T(m.reg.size * 4);
	Integer r((word)0, m.reg.size);
	unsigned k = AlmostInverse(r.reg, T, reg, reg.size, m.reg, m.reg.size);
	DivideByPower2Mod(r.reg, r.reg, k, m.reg, m.reg.size);
	return r;
}

word Integer::InverseMod(const word mod) const
{
    word g0 = mod, g1 = *this % mod;
    word v0 = 0, v1 = 1;
    word y;

	while (g1)
    {
		if (g1 == 1)
			return v1;
        y = g0 / g1;
        g0 = g0 % g1;
        v0 += y * v1;

		if (!g0)
			break;
		if (g0 == 1)
			return mod-v0;
		y = g1 / g0;
		g1 = g1 % g0;
		v1 += y * v0;
    }
	return 0;
}

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

Integer ModularArithmetic::Add(const Integer &a, const Integer &b) const
{
	if (a.reg.size==modulus.reg.size && b.reg.size==modulus.reg.size)
	{
		if (CryptoPP::Add(result.reg.ptr, a.reg, b.reg, a.reg.size)
			|| Compare(result.reg, modulus.reg, a.reg.size) >= 0)
		{
			CryptoPP::Subtract(result.reg.ptr, result.reg, modulus.reg, a.reg.size);
		}
		return result;
	}
	else
		{Integer r=a+b; if (r>=modulus) r-=modulus; return r;}
}

Integer& ModularArithmetic::Accumulate(Integer &a, const Integer &b) const
{
	if (a.reg.size==modulus.reg.size && b.reg.size==modulus.reg.size)
	{
		if (CryptoPP::Add(a.reg, a.reg, b.reg, a.reg.size)
			|| Compare(a.reg, modulus.reg, a.reg.size) >= 0)
		{
			CryptoPP::Subtract(a.reg, a.reg, modulus.reg, a.reg.size);
		}
	}
	else
		{a+=b; if (a>=modulus) a-=modulus;}

	return a;
}

Integer ModularArithmetic::Subtract(const Integer &a, const Integer &b) const
{
	if (a.reg.size==modulus.reg.size && b.reg.size==modulus.reg.size)
	{
		if (CryptoPP::Subtract(result.reg.ptr, a.reg, b.reg, a.reg.size))
			CryptoPP::Add(result.reg.ptr, result.reg, modulus.reg, a.reg.size);
		return result;
	}
	else
		return Add(a, Inverse(b));
}

Integer& ModularArithmetic::Reduce(Integer &a, const Integer &b) const
{
	if (a.reg.size==modulus.reg.size && b.reg.size==modulus.reg.size)
	{
		if (CryptoPP::Subtract(a.reg, a.reg, b.reg, a.reg.size))
			CryptoPP::Add(a.reg, a.reg, modulus.reg, a.reg.size);
	}
	else
		Accumulate(a, Inverse(b));

	return a;
}

Integer ModularArithmetic::Inverse(const Integer &a) const
{
	if (!a)
		return a;

	CopyWords(result.reg.ptr, modulus.reg, modulus.reg.size);
	if (CryptoPP::Subtract(result.reg.ptr, result.reg, a.reg, a.reg.size))
		Decrement(result.reg.ptr+a.reg.size, 1, modulus.reg.size-a.reg.size);

	return result;
}

Integer ModularArithmetic::MultiplicativeInverse(const Integer &a) const
{
	return a.InverseMod(modulus);
}

Integer ModularArithmetic::Exponentiate(const Integer &a, const Integer &e) const
{
	if (modulus.IsOdd())
	{
		MontgomeryRepresentation dr(modulus);
		return dr.ConvertOut(dr.Exponentiate(dr.ConvertIn(a), e));
	}
	else
		return AbstractRing<Integer>::Exponentiate(a, e);
}

Integer ModularArithmetic::CascadeExponentiate(const Integer &x, const Integer &e1, const Integer &y, const Integer &e2) const
{
	if (modulus.IsOdd())
	{
		MontgomeryRepresentation dr(modulus);
		return dr.ConvertOut(dr.CascadeExponentiate(dr.ConvertIn(x), e1, dr.ConvertIn(y), e2));
	}
	else
		return AbstractRing<Integer>::CascadeExponentiate(x, e1, y, e2);
}

MontgomeryRepresentation::MontgomeryRepresentation(const Integer &m)	// modulus must be odd
	: ModularArithmetic(m),
	  u((word)0, modulus.reg.size),
	  workspace(5*modulus.reg.size)
{
	assert(modulus.IsOdd());
	RecursiveInverseModPower2(u.reg, workspace, modulus.reg, modulus.reg.size);
}

Integer MontgomeryRepresentation::Multiply(const Integer &a, const Integer &b) const
{
	word *const T = workspace.ptr;
	word *const R = result.reg.ptr;
	const unsigned int N = modulus.reg.size;
	assert(a.reg.size<=N && b.reg.size<=N);

	AsymmetricMultiply(T, T+2*N, a.reg, a.reg.size, b.reg, b.reg.size);
	SetWords(T+a.reg.size+b.reg.size, 0, 2*N-a.reg.size-b.reg.size);
	MontgomeryReduce(R, T+2*N, T, modulus.reg, u.reg, N);
	return result;
}

Integer MontgomeryRepresentation::Square(const Integer &a) const
{
	word *const T = workspace.ptr;
	word *const R = result.reg.ptr;
	const unsigned int N = modulus.reg.size;
	assert(a.reg.size<=N);

	RecursiveSquare(T, T+2*N, a.reg, a.reg.size);
	SetWords(T+2*a.reg.size, 0, 2*N-2*a.reg.size);
	MontgomeryReduce(R, T+2*N, T, modulus.reg, u.reg, N);
	return result;
}

Integer MontgomeryRepresentation::ConvertOut(const Integer &a) const
{
	word *const T = workspace.ptr;
	word *const R = result.reg.ptr;
	const unsigned int N = modulus.reg.size;
	assert(a.reg.size<=N);

	CopyWords(T, a.reg, a.reg.size);
	SetWords(T+a.reg.size, 0, 2*N-a.reg.size);
	MontgomeryReduce(R, T+2*N, T, modulus.reg, u.reg, N);
	return result;
}

Integer MontgomeryRepresentation::MultiplicativeInverse(const Integer &a) const
{
//    return (EuclideanMultiplicativeInverse(a, modulus)<<(2*WORD_BITS*modulus.reg.size))%modulus;
	word *const T = workspace.ptr;
	word *const R = result.reg.ptr;
	const unsigned int N = modulus.reg.size;
	assert(a.reg.size<=N);

	CopyWords(T, a.reg, a.reg.size);
	SetWords(T+a.reg.size, 0, 2*N-a.reg.size);
	MontgomeryReduce(R, T+2*N, T, modulus.reg, u.reg, N);
	unsigned k = AlmostInverse(R, T, R, N, modulus.reg, N);

//	cout << "k=" << k << " N*32=" << 32*N << endl;

	if (k>N*WORD_BITS)
		DivideByPower2Mod(R, R, k-N*WORD_BITS, modulus.reg, N);
	else
		MultiplyByPower2Mod(R, R, N*WORD_BITS-k, modulus.reg, N);

	return result;
}

HalfMontgomeryRepresentation::HalfMontgomeryRepresentation(const Integer &m)	// modulus must be odd
	: ModularArithmetic(m),
	  v((modulus.reg.CleanGrow(4), Integer::Power2(3*WORD_BITS*modulus.reg.size/2)%modulus)),
	  u((word)0, modulus.reg.size/2),
	  workspace(4*modulus.reg.size)
{
	assert(modulus.IsOdd());

	result.reg.Grow(4);
	v.reg.CleanGrow(modulus.reg.size);
	RecursiveInverseModPower2(u.reg, workspace, modulus.reg, modulus.reg.size/2);
}

Integer HalfMontgomeryRepresentation::Multiply(const Integer &a, const Integer &b) const
{
	word *const T = workspace.ptr;
	word *const R = result.reg.ptr;
	const unsigned int N = modulus.reg.size;
	assert(a.reg.size<=N && b.reg.size<=N);

	AsymmetricMultiply(T, T+2*N, a.reg, a.reg.size, b.reg, b.reg.size);
	SetWords(T+a.reg.size+b.reg.size, 0, 2*N-a.reg.size-b.reg.size);
	HalfMontgomeryReduce(R, T+2*N, T, modulus.reg, u.reg, v.reg, N);
	return result;
}

Integer HalfMontgomeryRepresentation::Square(const Integer &a) const
{
	word *const T = workspace.ptr;
	const unsigned int N = modulus.reg.size;
	word *const R = result.reg.ptr;
	assert(a.reg.size<=N);

	RecursiveSquare(T, T+2*N, a.reg, a.reg.size);
	SetWords(T+2*a.reg.size, 0, 2*N-2*a.reg.size);
	HalfMontgomeryReduce(R, T+2*N, T, modulus.reg, u.reg, v.reg, N);
	return result;
}

Integer HalfMontgomeryRepresentation::ConvertOut(const Integer &a) const
{
	word *const T = workspace.ptr;
	word *const R = result.reg.ptr;
	const unsigned int N = modulus.reg.size;
	assert(a.reg.size<=N);

	CopyWords(T, a.reg, a.reg.size);
	SetWords(T+a.reg.size, 0, 2*N-a.reg.size);
	HalfMontgomeryReduce(R, T+2*N, T, modulus.reg, u.reg, v.reg, N);
	return result;
}

Integer HalfMontgomeryRepresentation::MultiplicativeInverse(const Integer &a) const
{
	return (a.InverseMod(modulus)<<(WORD_BITS*modulus.reg.size))%modulus;
}

NAMESPACE_END