/*
Two's-complement variable-sized bigInt number class by M Phillips - 2008.
Thanks also to Zero Soma Valintine for several bug fixes

This code is provided as is with no warranties or guarantees of
any kind.

Please send an email to M Phillips (mbp2@i4free.co.nz)
  - if you use this file in a released product, or
  - if you find any bugs, or
  - if you have any suggestions
*/

#include <assert.h>
#include <sstream>
#include <utility>
#include <algorithm>
#include <cmath>
#if defined(_MSC_VER) && _MSC_VER >= 1500
#include <intrin.h>
#endif
#include "varbigint.h"

#define BASE_MAX (std::numeric_limits<varbigint::BIGINT_BASE>::max)()

/* 'resize' resizes the dynamic array to the length 'newLen'. */
varbigint& varbigint::resize(int newLen) {
	if (newLen <= allocLen && newLen * 2 >= allocLen) {
		const BIGINT_BASE fill = negative ? BASE_MAX : 0U;
		for (int i = len; i < newLen; ++i)
			data[i] = fill;
		len = newLen;
	} else if (newLen != len) {
		varbigint temp(*this, newLen);
		swap(temp);
	}
	return *this;
}

/* 'trim' fixes up the length. */
varbigint& varbigint::trim() {
	const BIGINT_BASE fill = negative ? BASE_MAX : 0U;
	while (len > 0 && data[len-1] == fill)
		--len;
	if (len * 2 < allocLen)
		return resize(len);
	return *this;
}

/* 'verify' makes sure it is trimmed correctly. */
bool varbigint::verify() {
	return (len <= allocLen) && ((len == 0) || (data != NULL && data[len-1] != (negative ? BASE_MAX : 0U)));
}

varbigint::varbigint(const varbigint &other) : info(other.info), allocLen(other.len),
	data(other.len > 0 ? new BIGINT_BASE[other.len] : NULL) {
	for (int i=0; i < len; ++i)
		data[i] = other.data[i];
}

varbigint::varbigint(const varbigint &other, int len) : len(len), negative(other.negative),
    allocLen(len), data(len > 0 ? new BIGINT_BASE[len] : NULL) {
	if (other.len < len) {
		int i;
		for (i = 0; i < other.len; ++i)
			data[i] = other.data[i];
		const BIGINT_BASE fill = negative ? BASE_MAX : 0U;
		for (; i < len; ++i)
			data[i] = fill;
	} else {
		for (int i=0; i < len; ++i)
			data[i] = other.data[i];
	}
}

varbigint::varbigint(const std::string &s) : info(0), allocLen(0), data(NULL) {
	std::stringstream ss;
	ss << s.c_str();
	ss >> *this;
}

varbigint::varbigint(const std::wstring &ws) : info(0), allocLen(0), data(NULL) {
	std::stringstream wss;
	wss << ws.c_str();
	wss >> *this;
}

varbigint::varbigint(float f) : info(0), allocLen(0), data(NULL) {
	internalInitialiseFromFloating(f);
}

varbigint::varbigint(double d) : info(0), allocLen(0), data(NULL) {
	internalInitialiseFromFloating(d);
}

varbigint::varbigint(long double ld) : info(0), allocLen(0), data(NULL) {
	internalInitialiseFromFloating(ld);
}

varbigint::operator float() const {
	return internalExtractFromFloating<float>();
}

varbigint::operator double() const {
	return internalExtractFromFloating<double>();
}

varbigint::operator long double() const {
	return internalExtractFromFloating<long double>();
}

#if defined(_MSC_VER) && _MSC_VER >= 1600
varbigint::varbigint(varbigint &&other) :
	info(other.info), allocLen(other.allocLen), data(other.data) {
		other.data = NULL;
}

varbigint& varbigint::operator = (varbigint &&rhs) {
	using std::swap;
	swap(data, rhs.data);
	info = rhs.info;
	allocLen = rhs.allocLen;
	return *this;
}
#endif

/* The swap member function is more efficient than using
 three assignment statements because it doesn't require
 any dynamic memory allocation.*/
void varbigint::swap(varbigint &other) {
	using std::swap;
	swap(info, other.info);
	swap(data, other.data);
	swap(allocLen, other.allocLen);
}

varbigint& varbigint::operator = (const varbigint &rhs) {
	varbigint temp(rhs);
	swap(temp);
	return *this;
}

varbigint::operator std::string() const {
	std::stringstream ss;
	ss << *this;
	return ss.str();
}

varbigint::operator std::wstring() const {
	std::wstringstream wss;
	wss << *this;
	return wss.str();
}

std::ostream& operator << (std::ostream &os, const varbigint &val) {
	return varbigint::output(os, val);
}

std::wostream& operator << (std::wostream &wos, const varbigint &val) {
	return varbigint::output(wos, val);
}

std::istream& operator >> (std::istream &is, varbigint &val) {
	return varbigint::input(is, val);
}

std::wistream& operator >> (std::wistream &wis, varbigint &val) {
	return varbigint::input(wis, val);
}

void varbigint::bitClear(int b) {
	if (negative) {
		int requiredlen = 1 + b/BASE_BITS;
		if (requiredlen > len)
			resize(requiredlen);
	}
	if (b / BASE_BITS < len)
		data[b / BASE_BITS] &= ~(1<<(b % BASE_BITS));
	trim();
}

void varbigint::bitSet(int b) {
	if (!negative) {
		int requiredlen = 1 + b/BASE_BITS;
		if (requiredlen > len)
			resize(requiredlen);
	}
	if (b / BASE_BITS < len)
		data[b / BASE_BITS] |= static_cast<varbigint::BIGINT_BASE>(1)<<(b % BASE_BITS);
	trim();
}

void varbigint::bitToggle(int b) {
	int requiredlen = 1 + b/BASE_BITS;
	if (requiredlen > len)
		resize(requiredlen);
	data[b / BASE_BITS] ^= static_cast<varbigint::BIGINT_BASE>(1)<<(b % BASE_BITS);
	trim();
}

bool varbigint::bitTest(int b) const {
	if (1 + b/BASE_BITS > len)
		return negative;
	return (data[b / BASE_BITS] & (static_cast<varbigint::BIGINT_BASE>(1)<<(b % BASE_BITS))) != 0;
}

// 'lowBitSet' is true if the number held is odd
bool varbigint::lowBitSet() const {
	if (len == 0)
		return negative;
	return (data[0] & 1U) != 0;
}

// Special version of division that's fast but only works with small divisors
// 'a' must be positive. quotient may alias a
void DivMod_short(const varbigint &a, varbigint::BIGINT_SHORT b, varbigint &quotient, varbigint::BIGINT_SHORT *remainder) {
	if (a.negative)
		throw std::out_of_range("Unexpected negative number");
	if (!b) {
		assert(false);
		b = 1 / b;  // Force a divide by zero exception. (b == 0)
		return;
	}
	const varbigint::BIGINT_SHORT *u = (const varbigint::BIGINT_SHORT*)a.data;
	varbigint newQuotient(varbigint(), quotient.len);
	varbigint::BIGINT_SHORT *q = (varbigint::BIGINT_SHORT*)newQuotient.data;
	varbigint::BIGINT_BASE k = 0;
	int aLen = a.len << 1;
#ifdef BIGINT_USE_BIG_ENDIAN
	for (int j = 0; j < aLen; ++j)
#else
	for (int j = aLen - 1; j >= 0; --j)
#endif
	{
		varbigint::BIGINT_BASE fullWord = (k<<varbigint::HALF_BASE_BITS) + u[j];
		q[j] = static_cast<varbigint::BIGINT_SHORT>(fullWord / b);
		k = fullWord - q[j] * b;
	}
	quotient.trim();
	quotient.swap(newQuotient);
	if (remainder != NULL)
		*remainder = static_cast<varbigint::BIGINT_SHORT>(k);
}

/* 'DivMod' takes advantage of rules such as -a/-b equals a/b or
 a/-b equals -(a/b) etc, to remove the need to deal with negatives
 during the actual divide. */
// remainder may alias a or b
// but a and b must be positive
void DivMod(const varbigint &a, const varbigint &b, varbigint &quotient, varbigint *remainder) {
	if (a.negative) {
		if (b.negative)
			DivMod_Unsigned(-a, -b, quotient, remainder);
		else {
			DivMod_Unsigned(-a, b, quotient, remainder);
			quotient = -quotient;
		}
		if (remainder != NULL) remainder->makeNegative();
	} else {
		if (b.negative) {
			DivMod_Unsigned(a, -b, quotient, remainder);
			quotient = -quotient;
		} else
			DivMod_Unsigned(a, b, quotient, remainder);
	}
}

// quotient may alias a or b, remainder may also alias a or b
// but a and b must be positive
void DivMod_Unsigned(const varbigint &a, const varbigint &b, varbigint &quotient, varbigint *remainder) {
	assert(!a.negative && !b.negative);
	int shiftcnt = 0;

	// Check for attempt to divide by zero
	if (!b) {
		shiftcnt = 1 / shiftcnt;  // Force a divide by zero exception. (shiftcnt=0)
	} else {
		varbigint a2 = a, b2 = b, c;
		// Left shift B until it is the same magnitude as A
		shiftcnt = findHighestBitSet(a2) - findHighestBitSet(b2);
		if (shiftcnt > 0)
			b2 <<= shiftcnt;

		if (b2 > a2) {		// If B is greater than A, right shift B
			b2 >>= 1;
			--shiftcnt;
		}
		while (shiftcnt >= 0) {
			int result = compare(a2, b2);
			if (result >= 0) {	// If B is not greater than A, then the bit is a 1
				a2 -= b2;		// Subtract B from A
				c.bitSet(shiftcnt);
				if (!result)	// Bail early if the remainder reaches zero
					break;
			}
			b2 >>= 1;		// Right shift B
			--shiftcnt;
		}

		quotient.swap(c);
		if (remainder != NULL)
			a2.swap(*remainder);
	}
}

varbigint& varbigint::operator &=(const varbigint &rhs) {
	int theLen = std::min(len, rhs.len);
	if ((rhs.len <= len && rhs.negative) || (len <= rhs.len && !negative)) {
		for (int i = 0; i < theLen; ++i)
			data[i] &= rhs.data[i];
		trim();
	} else {
		varbigint temp(rhs);
		for (int i = 0; i < theLen; ++i)
			temp.data[i] &= data[i];
		temp.trim();
		swap(temp);
	}
	return *this;
}

varbigint& varbigint::operator |=(const varbigint &rhs) {
	int theLen = std::min(len, rhs.len);
	if ((rhs.len <= len && !rhs.negative) || (len <= rhs.len && negative)) {
		for (int i = 0; i < theLen; ++i)
			data[i] |= rhs.data[i];
		trim();
	} else {
		varbigint temp(rhs);
		for (int i = 0; i < theLen; ++i)
			temp.data[i] |= data[i];
		temp.trim();
		swap(temp);
	}
	return *this;
}

varbigint& varbigint::operator ^=(const varbigint &rhs) {
	if (rhs.len <= len) {
		for (int i = 0; i < rhs.len; ++i)
			data[i] ^= rhs.data[i];
		negative ^= rhs.negative;
		trim();
	} else {
		varbigint temp(rhs);
		for (int i = 0; i < len; ++i)
			temp.data[i] ^= data[i];
		temp.negative ^= negative;
		temp.trim();
		swap(temp);
	}
	return *this;
}

/* 'operator >>=' works out exactly where the bits should end up
 and puts them straight there. This is faster than rotating by
 one bit successively. Lower bits get thrown away, and higher
 bits get filled with zero or one depending on the sign. */
varbigint& varbigint::operator>>=(int shift) {
	int source = shift / BASE_BITS;
	int remaindershift = shift & (BASE_BITS-1);
	int othershift = BASE_BITS - remaindershift;
	const BIGINT_BASE fill = negative ? BASE_MAX : 0U;

	for (int i = 0; i < len; ++i) {
		if (source < len) {
			data[i] = data[source] >> remaindershift;
			if (++source < len && othershift < BASE_BITS)
				data[i] |= data[source] << othershift;
			else
				data[i] |= fill << othershift;
		} else {
			data[i] = fill;
		}
	}
	return trim();
}

/* 'operator <<=' works just like the >>= operator, except that
 lower bits get zero filled, and higher bits just move over. */
varbigint& varbigint::operator<<=(int shift) {
	resize(len + (shift+BASE_BITS-1)/BASE_BITS);

	int source = len-1 - shift / BASE_BITS;
	int remaindershift = shift & (BASE_BITS-1);
	int othershift = BASE_BITS - remaindershift;

	for (int i = len-1; i >= 0; --i) {
		if (source >= 0) {
			data[i] = data[source] << remaindershift;
			if (--source >= 0 && othershift < BASE_BITS)
				data[i] |= data[source] >> othershift;
		} else {
			data[i] = 0U;
		}
	}
	return trim();
}

varbigint& varbigint::operator +=(const varbigint &rhs) {
	if (negative) {
		varbigint a(-*this);
		if (rhs.negative)
			a += -rhs;
		else
			a -= rhs;
		return *this = -a;
	} else if (rhs.negative) {
		return *this -= -rhs;
	}
	varbigint a(*this);
	varbigint b(rhs);
	if (a.len < b.len) {
		a.swap(b);
	}
	BIGINT_BASE carry = 0U;
	for (int i = 0; i < b.len; ++i) {
		a.data[i] += b.data[i] + carry;
		if (!carry) {
			carry = (a.data[i] <  b.data[i]);
		} else {
			carry = (a.data[i] <= b.data[i]);
		}
	}
	if (carry) {
		int i = b.len;
		do {
			if (i >= a.len)
				a.resize(i+1);
		} while (!++a.data[i++]);
	}
	a.trim();
	swap(a);
	return *this;
}

varbigint& varbigint::operator -=(const varbigint &rhs) {
	if (negative) {
		varbigint a(-*this);
		if (rhs.negative)
			a -= -rhs;
		else
			a += rhs;
		return *this = -a;
	} else if (rhs.negative) {
		return *this += -rhs;
	}
	varbigint a(*this);
	varbigint b(rhs);
	bool swapped = false;
	if (a < b) {
		a.swap(b);
		swapped = true;
	}
	BIGINT_BASE borrow = 0U, prevdigit;
	for (int i = 0; i < b.len; ++i) {
		prevdigit = a.data[i];
		a.data[i] -= b.data[i] + borrow;
		if (!borrow) {
			borrow = (prevdigit <  b.data[i]);
		} else {
			borrow = (prevdigit <= b.data[i]);
		}
	}
	if (borrow) {
		int i = b.len;
		do {
			if (i >= a.len)
				a.resize(i+1);
		} while (--a.data[i++] == BASE_MAX);
	}
	if (swapped)
		a.makeNegative();
	a.trim();
	swap(a);
	return *this;
}

/* 'Multiply' takes advantage of rules such as -a*-b equals a*b or
 a*-b equals -(a*b) etc, to remove the need to deal with negatives
 during the actual multiply. */
void varbigint::Multiply(const varbigint &a, const varbigint &b, varbigint &result) {
	if (a.negative) {
		if (b.negative)
			varbigint::Multiply_Unsigned(-a, -b, result);
		else {
			varbigint::Multiply_Unsigned(-a, b, result);
			result.makeNegative();
		}
	} else {
		if (b.negative) {
			varbigint::Multiply_Unsigned(a, -b, result);
			result.makeNegative();
		} else
			varbigint::Multiply_Unsigned(a, b, result);
	}
}

void varbigint::Multiply_Unsigned(const varbigint &a, const varbigint &b, varbigint &result) {
	assert(!a.negative && !b.negative);
	int aLen = a.len, bLen = b.len;
	if (aLen > 175 && bLen > 175) {
		// Karatsuba steps
		int shift = (aLen + bLen) * (BASE_BITS >> 2);
		varbigint lowMask(0);
		lowMask.bitSet(shift);
		--lowMask;
		varbigint x0 = a & lowMask;
		varbigint x1 = a >> shift;
		varbigint y0 = b & lowMask;
		varbigint y1 = b >> shift;
		varbigint z2 = x1 * y1;
		varbigint z0 = x0 * y0;
		varbigint z1 = (x1 + x0) * (y1 + y0) - z2 - z0;
		result = (z2 << (shift<<1)) + (z1 << shift) + z0;
		return;
	}
	varbigint res(varbigint(), aLen + bLen);
	aLen <<= 1; bLen <<= 1;
	BIGINT_SHORT *w = (BIGINT_SHORT*)res.data;	// w, u, v reference the
	const BIGINT_SHORT *u = (BIGINT_SHORT*)a.data;	// fullword arrays as
	const BIGINT_SHORT *v = (BIGINT_SHORT*)b.data;	// halfwords.
#ifdef BIGINT_USE_BIG_ENDIAN
	for (int j = bLen-1; j >= 0; --j) {
		BIGINT_BASE k = 0, vj = v[j];
		for (int i = aLen-1; i >= 0; --i) {
			const int ij = i + j + 1;
			k += u[i]*vj + w[ij];
			w[ij] = static_cast<varbigint::BIGINT_SHORT>(k);
			k >>= HALF_BASE_BITS;
		}
		w[j] = static_cast<varbigint::BIGINT_SHORT>(k);
	}
#else
	for (int j = 0; j < bLen; ++j) {
		BIGINT_BASE k = 0, vj = v[j];
		for (int i = 0; i < aLen; ++i) {
			const int ij = i + j;
			k += u[i]*vj + w[ij];
			w[ij] = static_cast<varbigint::BIGINT_SHORT>(k);
			k >>= HALF_BASE_BITS;
		}
		w[aLen+j] = static_cast<varbigint::BIGINT_SHORT>(k);
	}
#endif
	res.trim();
	result.swap(res);
}

varbigint& varbigint::Mul_short(const varbigint::BIGINT_SHORT rhs) {
	varbigint result(varbigint(), len + 1);
	int aLen = len << 1;
	varbigint::BIGINT_SHORT *w = (varbigint::BIGINT_SHORT*)result.data;		// w, u reference the
	const varbigint::BIGINT_SHORT *u = (const varbigint::BIGINT_SHORT*)data;	// fullword arrays as halfwords
	varbigint::BIGINT_BASE k = 0;
#ifdef BIGINT_USE_BIG_ENDIAN
	for (int i = aLen-1; i >= 0; --i) {
		k += u[i]*rhs + w[i];
		w[i+2] = static_cast<varbigint::BIGINT_SHORT>(k);
		k >>= varbigint::HALF_BASE_BITS;
	}
	result.data[0] = k;
#else
	for (int i = 0; i < aLen; ++i) {
		k += u[i]*rhs + w[i];
		w[i] = static_cast<varbigint::BIGINT_SHORT>(k);
		k >>= varbigint::HALF_BASE_BITS;
	}
	result.data[len] = k;
#endif
	result.trim();
	swap(result);
	return *this;
}

varbigint& varbigint::Mul_int(const varbigint::BIGINT_BASE rhs) {
	varbigint result(varbigint(), len + 1);
	int aLen = len << 1;
	BIGINT_SHORT *w = (BIGINT_SHORT*)result.data;	// w, u, v reference the
	const BIGINT_SHORT *u = (BIGINT_SHORT*)data;	// fullword arrays as
	const BIGINT_SHORT *v = (BIGINT_SHORT*)&rhs;	// halfwords.
#ifdef BIGINT_USE_BIG_ENDIAN
	for (int j = 1; j >= 0; --j) {
		BIGINT_BASE k = 0, vj = v[j];
		for (int i = aLen-1; i >= 0; --i) {
			const int ij = i + j + 1;
			k += u[i]*vj + w[ij];
			w[ij] = static_cast<varbigint::BIGINT_SHORT>(k);
			k >>= HALF_BASE_BITS;
		}
		w[j] = static_cast<varbigint::BIGINT_SHORT>(k);
	}
#else
	for (int j = 0; j < 2; ++j) {
		BIGINT_BASE k = 0, vj = v[j];
		for (int i = 0; i < aLen; ++i) {
			const int ij = i + j;
			k += u[i]*vj + w[ij];
			w[ij] = static_cast<varbigint::BIGINT_SHORT>(k);
			k >>= HALF_BASE_BITS;
		}
		w[aLen+j] = static_cast<varbigint::BIGINT_SHORT>(k);
	}
#endif
	result.trim();
	swap(result);
	return *this;
}

varbigint& varbigint::operator *=(const varbigint &rhs) {
	Multiply(*this, rhs, *this);
	return *this;
}

varbigint& varbigint::operator /=(const varbigint &rhs) {
	DivMod(*this, rhs, *this, NULL);
	return *this;
}

varbigint& varbigint::operator %=(const varbigint &rhs) {
	varbigint quotient;
	DivMod(*this, rhs, quotient, this);
	return *this;
}

int compare(const varbigint &a, const varbigint &b) {
	if (a.negative != b.negative)
		return a.negative ? -1 : 1;
	if (a.len != b.len)
		return (a.len < b.len) != static_cast<bool>(a.negative) ? -1 : 1;
	for (int i = a.len-1; i >= 0; --i) {
		if (a.data[i] < b.data[i])
			return -1;
		if (a.data[i] > b.data[i])
			return 1;
	}
	return 0;
}

/* 'operator ==' must check three things for equality:
 the sign, the length, and the array data (in that order).
*/
bool operator ==(const varbigint &a, const varbigint &b) {
	if (a.negative != b.negative)
		return false;
	if (a.len != b.len)
		return false;
	for (int i = 0; i < a.len; ++i)
		if (a.data[i] != b.data[i])
			return false;
	return true;
}

bool operator !=(const varbigint &a, const varbigint &b) {
	return !(a == b);
}

bool operator < (const varbigint &a, const varbigint &b) {
	if (a.negative != b.negative)
		return a.negative;
	if (a.len != b.len)
		return (a.len < b.len) != static_cast<bool>(a.negative);
	for (int i = a.len-1; i >= 0; --i) {
		if (a.data[i] < b.data[i])
			return true;
		if (a.data[i] > b.data[i])
			return false;
	}
	return false;
}

bool operator > (const varbigint &a, const varbigint &b) {
	return (b < a);
}

bool operator <= (const varbigint &a, const varbigint &b) {
	return !(b < a);
}

bool operator >= (const varbigint &a, const varbigint &b) {
	return !(a < b);
}

const varbigint varbigint::operator>> (int shift) const {
	varbigint result = *this;
	result >>= shift;
	return result;
}

const varbigint varbigint::operator<< (int shift) const {
	varbigint result = *this;
	result <<= shift;
	return result;
}

const varbigint operator & (varbigint a, const varbigint &b) {
	return a &= b;
}

const varbigint operator | (varbigint a, const varbigint &b) {
	return a |= b;
}

const varbigint operator ^ (varbigint a, const varbigint &b) {
	return a ^= b;
}

const varbigint operator + (varbigint a, const varbigint &b) {
	return a += b;
}

const varbigint operator - (varbigint a, const varbigint &b) {
	return a -= b;
}

const varbigint operator * (const varbigint &a, const varbigint &b) {
	varbigint result;
	varbigint::Multiply(a, b, result);
	return result;
}

const varbigint operator / (const varbigint &a, const varbigint &b) {
	varbigint quotient;
	DivMod(a, b, quotient, NULL);
	return quotient;
}

const varbigint operator % (const varbigint &a, const varbigint &b) {
	varbigint quotient, result;
	DivMod(a, b, quotient, &result);
	return result;
}

varbigint& varbigint::operator++ () {  // Pre Increment operator -- faster than add
	// We must first check for all bits being set, which is a special case
	// because we either need to change the sign, or we need to extend the array
	// because we've reached a critical power of two that requires more storage.
	bool result = true;
	for (int i = 0; i < len; ++i) {
		if (data[i] != BASE_MAX) {
			result = false;
			break;
		}
	}
	if (result) {
		if (negative) {
			negative = false;
			return resize(0);
		}
		resize(len+1);
	}
	// Now we perform an add of one and deal with carries
	int i = 0;
	while (!++data[i++]) {}
	return trim();
}

const varbigint varbigint::operator++ (int) {  // Post Increment operator -- faster than add
	varbigint result = *this;
	++*this;
	return result;
}

varbigint& varbigint::operator-- () {  // Pre Decrement operator -- faster than subtract
	// We must first check for all bits being clear, which is a special case
	// because we either need to change the sign, or we need to extend the array
	// because we've reached a critical negative power of two that requires more
	// storage.
	bool result = true;
	for (int i = 0; i < len; ++i) {
		if (data[i] != 0) {
			result = false;
			break;
		}
	}
	if (result) {
		if (!negative) {
			negative = true;
			return resize(0);
		}
		resize(len+1);
	}
	// Now we perform a subtract of one and deal with borrows
	int i = 0;
	while (--data[i++] == BASE_MAX) {}
	return trim();
}

const varbigint varbigint::operator-- (int) {  // Post Decrement operator -- faster than subtract
	varbigint result = *this;
	--*this;
	return result;
}

bool varbigint::operator ! () const {	//For comparison against zero
	if (negative)
		return false;
	for (int i = 0; i < len; ++i)
		if (data[i])
			return false;
	return true;
}

/* 'operator ~' returns a value with all the bits flipped AND the
 sign flipped. */
const varbigint varbigint::operator ~ () const {
	varbigint result;
	result.resize(len);
	for (int i = 0; i < len; ++i)
		result.data[i] = ~data[i];
	result.negative = !negative;
	return result;
}

const varbigint& varbigint::operator + () const {  // Unary positive
	return *this;
}

/* 'operator -' performs standard two's-complement negation.
 Either flips bits then adds one, or subtracts one then flips bits.
*/
const varbigint varbigint::operator - () const {  // Unary Negative
	varbigint result;
	if (!negative) {
		result.resize(len);
		result.negative = true;
		for (int i = 0; i < len; ++i)
			result.data[i] = ~data[i];
		++result;
	} else {
		result = *this;
		--result;
		for (int i = 0; i < result.len; ++i)
			result.data[i] = ~result.data[i];
		result.negative = !negative;
		result.trim();
	}
	return result;
}

void varbigint::makeNegative() {
	if (!negative) {
		negative = true;
		for (int i = 0; i < len; ++i)
			data[i] = ~data[i];
		++*this;
	}
}

void varbigint::makePositive() {
	if (negative) {
		--*this;
		for (int i = 0; i < len; ++i)
			data[i] = ~data[i];
		negative = !negative;
		trim();
	}
}

/* 'sqrt' finds the square root of a number by the Newton-Raphson
method. */
const varbigint sqrt(const varbigint &x) {
	if (x.negative)
		return 0U;	// if negative just return zero
	int bit = findHighestBitSet(x);
	// "1<<bit" starts at the highest power of four <= the argument.
	bit -= (bit & 1);
	varbigint num = x, res = 0;

	while (bit >= 0) {
		varbigint temp = res;
		temp.bitSet(bit);
		res >>= 1;
		if (num >= temp) {
			num -= temp;
			res.bitSet(bit);
		}
		bit -= 2;
	}
	return res;
}

/* 'cbrt' finds the cube root of a number by the Newton-Raphson
method. */
const varbigint cbrt(const varbigint &x) {
	varbigint x1, y1, b, y2, y3;
	bool negative = false;
	if (x < 0) {
		negative = true;
		x1 = -x;
	} else
		x1 = x;
	int bit = (findHighestBitSet(x1) / 3) * 3;
	// "1<<bit" starts at the highest power of eight <= the argument.
	y1 = 0;
	y2 = 0;
	while (bit >= 0) {
		y2 <<= 2;
		y1 <<= 1;
		y3 = y2 + y1;
		b = ((y3<<1) + y3 + 1) << bit;
		if (x1 >= b) {
			x1 -= b;
			y2 += (y1<<1) + 1;
			++y1;
		}
		bit -= 3;
	}
	return negative ? -y1 : y1;
}

const varbigint abs(const varbigint &rhs) {
	return rhs.negative ? -rhs : rhs;
}

const varbigint factorial(const varbigint &rhs) {
	if (rhs.negative)
		return 0;
	varbigint result = 1;
	varbigint::BIGINT_BASE f = rhs;
	//If the number to take the factorial of is bigger than can fit into
	//an unsigned long then the there's no way in hell the result can be
	//stored in conventional PC memory anyway - TOO BIG!
	while (f > static_cast<varbigint::BIGINT_BASE>(varbigint::LOWER_HALF_MASK))
		result.Mul_int(f--);
	while (f > 0U)
		result.Mul_short((varbigint::BIGINT_SHORT)f--);
	return result;
}

int jacobi(varbigint m, varbigint n) {
	assert(!m.negative && !n.negative);
	assert(n.lowBitSet());
	int sign = 1;
	for (;;) {
		m %= n;
		if (!m)
			return 0;
		int lo = findLowestBitSet(m);
		if (lo > 0) {
			m >>= lo;
			if ((lo & 1) != 0 && ((1 << static_cast<int>(n&7)) & ((1<<3) | (1<<5))) != 0)
				sign = -sign;
		}
		if (m == 1)
			return sign;
		if ((3 == static_cast<int>(n&3)) && (3 == static_cast<int>(m&3)))
			sign = -sign;
		m.swap(n);
	}
}

static const unsigned int primesUnder1000[] = {
	  2,   3,   5,   7,  11,  13,  17,  19,  23,  29,
	 31,  37,  41,  43,  47,  53,  59,  61,  67,  71,
	 73,  79,  83,  89,  97, 101, 103, 107, 109, 113,
	127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
	179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
	233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
	283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
	353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
	419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
	467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
	547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
	607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
	661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
	739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
	811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
	877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
	947, 953, 967, 971, 977, 983, 991, 997,
};

bool Fermat(const varbigint &v, unsigned int maxa) {
	const unsigned int *end = primesUnder1000 + sizeof(primesUnder1000)/sizeof(primesUnder1000[0]);
	for (const unsigned int *it = primesUnder1000; it != end; ++it) {
		unsigned int a = *it;
		varbigint x = modpow(a, v, v);
		if (x != (int)a)
			return false;
		if (a >= maxa)
			break;
	}
	return true;
}

bool MillerRabin(const varbigint &v, unsigned int maxa) {
	varbigint vMinus1 = v - 1;
	unsigned int s = findLowestBitSet(vMinus1);
	varbigint d = vMinus1 >> (int)s;

	const unsigned int *end = primesUnder1000 + sizeof(primesUnder1000)/sizeof(primesUnder1000[0]);
	for (const unsigned int *it = primesUnder1000; it != end; ++it) {
		unsigned int a = *it;
		varbigint x = modpow(a, d, v);
		if (x != 1 && x != vMinus1) {
			unsigned int r = 1;
			varbigint dd = d;
			for (; r < s; r++) {
				dd <<= 1;
				x = modpow(a, dd, v);
				if (x == vMinus1)
					break;
			}
			if (r >= s)
				return false;
		}
		if (a >= maxa)
			break;
	}
	return true;
}

bool isPrime(const varbigint &v) {
	const unsigned int *end = primesUnder1000 + sizeof(primesUnder1000)/sizeof(primesUnder1000[0]);
	if (v <= 2 || !v.lowBitSet()) {
		return v == 2;
	} else if (v < 1000) {
		unsigned int uiv = v;
		return std::binary_search(primesUnder1000, end, uiv);
	} else if (v < 1000000) {
		unsigned int uiv = v;
		for (const unsigned int *it = primesUnder1000; it != end && ((*it) * (*it) <= uiv); ++it)
			if (!(uiv % *it))
				return false;
		return true;
	} else {
		for (const unsigned int *it = primesUnder1000; it != end; ++it)
			if (!(v % varbigint(*it)))
				return false;
		if (v < varbigint(1373653)) {
			return MillerRabin(v, 3);
		} else if (v < varbigint(25326001)) {
			return MillerRabin(v, 5);
		} else if (v < varbigint("3215031751")) {
			return MillerRabin(v, 7);
		} else if (v < varbigint("2152302898747")) {
			return MillerRabin(v, 11);
		} else if (v < varbigint("3474749660383")) {
			return MillerRabin(v, 13);
		} else if (v < varbigint("341550071728321")) {
			return MillerRabin(v, 17);
		} else if (v < varbigint("41234316135705689041")) {
			return MillerRabin(v, 23);
		} else if (v < varbigint("1553360566073143205541002401")) {
			return MillerRabin(v, 29);
		} else if (v < varbigint("56897193526942024370326972321")) {
			return MillerRabin(v, 31);
		} else {
		// Good luck if you end up hitting this case, it will take a long time.
		// Many further improvements are required here
			varbigint stop = sqrt(v);
			for (varbigint factor = 1009U; factor <= stop; factor+=2U) {
				if (!(v % factor)) {
					return false;
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
				}
			}
			return true;
		}
	}
}

/* 'gcd' is the Greatest Common Divisor algorithm. */
const varbigint gcd(const varbigint &a, const varbigint &b) {
	if (!a)
		return b;
	if (!b)
		return a;
	varbigint a2 = a, b2 = b;
	a2.makePositive();
	b2.makePositive();
	int aLowBitPos = findLowestBitSet(a2);
	int bLowBitPos = findLowestBitSet(b2);
	int shift = aLowBitPos < bLowBitPos ? aLowBitPos : bLowBitPos;
	a2 >>= aLowBitPos;
	b2 >>= bLowBitPos;
	do {
		while (!a2.lowBitSet())
			a2 >>= 1;
		if (b2 < a2) {
			a2 -= b2;
		} else {
			varbigint c = abs(b2 - a2);
			b2 = a2;
			a2 = c;
		}
		a2 >>= 1;
	} while (!!a2);
	return b2 << shift;
}

const varbigint lcm(const varbigint &a, const varbigint &b) {
	return (a / gcd(a, b)) * b;
}

const varbigint gcdext(varbigint a, varbigint b, varbigint &s, varbigint &t) {
	varbigint nexts = 0, nextt = 1;
	s = 1;
	t = 0;
	while (!!b) {
		varbigint quo, rem;
		DivMod(a, b, quo, &rem);
		b.swap(rem);
		a.swap(rem);
		varbigint s1 = s - quo*nexts;
		varbigint t1 = t - quo*nextt;
		nexts.swap(s1);
		s.swap(s1);
		nextt.swap(t1);
		t.swap(t1);
	}
	return a;
}

int findHighestBitSet(const varbigint &x) {
	int result = x.len*varbigint::BASE_BITS-1;
	for (int i = x.len-1; i >= 0; --i) {
		if (x.data[i] != 0) {
			varbigint::BIGINT_BASE temp = x.data[i];
			varbigint::BIGINT_BASE tempMax = BASE_MAX;
			int shift = varbigint::HALF_BASE_BITS;
			do {
				tempMax <<= shift;
				if ((temp & tempMax) == 0) {
					result -= shift;
					temp <<= shift;
				}
			} while ((shift>>=1) > 0);
			break;
		}
		result -= varbigint::BASE_BITS;
	}
	return result;
}

int findLowestBitSet(const varbigint &x) {
	int result = -1;
	for (int i = 0; i < x.len; ++i) {
		if (x.data[i] != 0) {
			++result;
			varbigint::BIGINT_BASE temp = x.data[i];
			varbigint::BIGINT_BASE tempMax = BASE_MAX;
			int shift = varbigint::HALF_BASE_BITS;
			do {
				tempMax >>= shift;
				if ((temp & tempMax) == 0) {
					result += shift;
					temp >>= shift;
				}
			} while ((shift>>=1) > 0);
			break;
		}
		result += varbigint::BASE_BITS;
	}
	return result;
}

bool isPow2(const varbigint &x) {
	bool found = false;
	for (int i = 0; i < x.len; ++i) {
		if (x.data[i] != 0) {
			if (found)
				return false;
			if ((x.data[i] & (x.data[i]-1)) != 0)
				return false;
			found = true;
		}
	}
	return true;
}

/* 'nextPow2' calculates the next highest number that is a
 power of two. */
const varbigint nextPow2(const varbigint &x) {
	int high = findHighestBitSet(x);
	if (high < 0)
		return 1;
	varbigint temp;
	temp.bitSet(high);
	return (temp < x) ? (temp << 1) : temp;
}

/* 'ceillog2' calculates the log base two of a number, rounded up. */
int ceillog2(const varbigint &x) {
	int high = findHighestBitSet(x);
	if (high < 0)
		return 0;
	varbigint temp;
	temp.bitSet(high);
	return (temp < x) ? (high + 1) : high;
}

template<typename T> struct varbigintPopCounter {
	static int popCount(T v) {
		int count = 0;
		while (v) {
			v &= v - 1; // clear the least significant bit set
			++count;
		}
		return count;
	}
};

#if defined(__GNUG__)
template<> struct varbigintPopCounter<unsigned int> {
	static int popCount(unsigned int v) { return __builtin_popcount(v); }
};
template<> struct varbigintPopCounter<unsigned long long> {
	static int popCount(unsigned long long v) { return static_cast<int>(__builtin_popcountll(v)); }
};
#elif defined(_MSC_VER) && _MSC_VER >= 1500 // >= VS2008
template<> struct varbigintPopCounter<unsigned int> {
	static int popCount(unsigned int v) { return __popcnt(v); }
};
template<> struct varbigintPopCounter<unsigned long long> {
#if _M_AMD64
	static int popCount(unsigned long long v) { return static_cast<int>(__popcnt64(v)); }
#else
	static int popCount(unsigned long long v) { return __popcnt((unsigned int)v) + __popcnt((unsigned int)(v>>32)); }
#endif
};
#endif

/* 'popcount' calculates how many bits are set in the number. */
int popcount(const varbigint &x) {
	if (x.negative)
		return 0;	// There are an infinite number of leading 1's - abort
	int count = 0;
	for (int i = 0; i < x.len; ++i)
		count += varbigintPopCounter<varbigint::BIGINT_BASE>::popCount(x.data[i]);
	return count;
}

/* 'pow' calculates a to the power of b using powers of two of 'a'*/
const varbigint pow(const varbigint &a, int b) {
	if (b < 0)
		return (varbigint)0;
	if (a == (varbigint)1U)
		return a;
	varbigint a2 = a;
	varbigint result = (varbigint)1U;
	while (b) {
		if (b & 1)
			result *= a2;
		b >>= 1;
		a2 *= a2;
	}
	return result;
}

const varbigint modpow(varbigint base, const varbigint &exp, const varbigint &mod)
{
	varbigint result = 1;
	int expBit = 0, lastExpBit = findHighestBitSet(exp);
	while (expBit <= lastExpBit) {
		if (exp.bitTest(expBit))
			result = (result * base) % mod;
		++expBit;
		base = (base * base) % mod;
	}
	return result;
}

#undef BASE_MAX