/***************************************************************************
 *   Copyright (C) 2005 by Jeff Ferr                                       *
 *   root@sat                                                              *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
#include "jmatrix.h"
#include "jmathlib.h"

namespace jmath {

BigInteger(int *val):
	jcommon::Object()
{
	jcommon::Object::SetClassName("jmath::BigInteger");

	if (val == NULL)
		throw new NumberFormatException("Zero length BigInteger");
	
	if (val[0] < 0) {
		mag = makePositive(val);
		signum = -1;
	} else {
		mag = trustedStripLeadingZeroInts(val);
		signum = (mag.length == 0 ? 0 : 1);
	}
}

BigInteger(int[] magnitude, int signum):
	jcommon::Object()
{
	jcommon::Object::SetClassName("jmath::BigInteger");
	
	this.signum = (magnitude.length==0 ? 0 : signum);
	this.mag = magnitude;
}

BigInteger(char[] magnitude, int signum):
	jcommon::Object()
{
	jcommon::Object::SetClassName("jmath::BigInteger");
	
	this.signum = (magnitude.length==0 ? 0 : signum);
	this.mag = stripLeadingZeroBytes(magnitude);
}

BigInteger(MutableBigInteger val, int sign):
	jcommon::Object()
{
	jcommon::Object::SetClassName("jmath::BigInteger");
	
	if (val.offset > 0 || val.value.length != val.intLen) {
		mag = new int[val.intLen];
		for(int i=0; i<val.intLen; i++)
			mag[i] = val.value[val.offset+i];
	} else {
		mag = val.value;
	}
	
	this.signum = (val.intLen == 0) ? 0 : sign;
}

BigInteger(long val):
	jcommon::Object()
{
	jcommon::Object::SetClassName("jmath::BigInteger");
	
	if (val < 0) {
		signum = -1;
		val = -val;
	} else {
		signum = 1;
	}
	
	int highWord = (int)(val >>> 32);
	if (highWord==0) {
		mag = new int[1];
		mag[0] = (int)val;
	} else {
		mag = new int[2];
		mag[0] = highWord;
		mag[1] = (int)val;
	}
}


BigInteger(int signum, int *magnitude):
	jcommon::Object()
{
	jcommon::Object::SetClassName("jmath::BigInteger");
	
	this.mag = stripLeadingZeroInts(magnitude);
	
	if (signum < -1 || signum > 1)
		throw(new NumberFormatException("Invalid signum value"));
	
	if (this.mag.length==0) {
		this.signum = 0;
	} else {
		if (signum == 0)
			throw(new NumberFormatException("signum-magnitude mismatch"));
		this.signum = signum;
	}
}


BigInteger(std::string val):
	jcommon::Object()
{
	jcommon::Object::SetClassName("jmath::BigInteger");
	
	this(val, 10);
}

BigInteger(int numBits, Random rnd):
	jcommon::Object()
{
	jcommon::Object::SetClassName("jmath::BigInteger");
	
	this(1, randomBits(numBits, rnd));
}

public BigInteger(int bitLength, int certainty, Random rnd):
	jcommon::Object()
{
	jcommon::Object::SetClassName("jmath::BigInteger");
	
	BigInteger prime;
	
	if (bitLength < 2)
		throw new ArithmeticException("bitLength < 2");
	// The cutoff of 95 was chosen empirically for best performance
	prime = (bitLength < 95 ? smallPrime(bitLength, certainty, rnd)
			: largePrime(bitLength, certainty, rnd));
	signum = 1;
	mag = prime.mag;
}


BigInteger(char *val):
	jcommon::Object()
{
	jcommon::Object::SetClassName("jmath::BigInteger");
	
	if (val == NULL)
		throw new NumberFormatException("Zero length BigInteger");
	
	if (val[0] < 0) {
		mag = makePositive(val);
		signum = -1;
	} else {
		mag = stripLeadingZeroBytes(val);
		signum = (mag.length == 0 ? 0 : 1);
	}
}

BigInteger(int signum, char *magnitude):
	jcommon::Object()
{
	jcommon::Object::SetClassName("jmath::BigInteger");
	
	this.mag = stripLeadingZeroBytes(magnitude);
	
	if (signum < -1 || signum > 1)
		throw(new NumberFormatException("Invalid signum value"));
	
	if (this.mag.length==0) {
		this.signum = 0;
	} else {
		if (signum == 0)
			throw(new NumberFormatException("signum-magnitude mismatch"));
		this.signum = signum;
	}
}

BigInteger(std::string val, int radix):
	jcommon::Object()
{
	jcommon::Object::SetClassName("jmath::BigInteger");
	
	int cursor = 0, numDigits;
	int len = val.length();
	
	if (radix < Character.MIN_RADIX || radix > Character.MAX_RADIX)
		throw new NumberFormatException("Radix out of range");
	if (val.length() == 0)
		throw new NumberFormatException("Zero length BigInteger");
	
	// Check for minus sign
	signum = 1;
	int index = val.lastIndexOf("-");
	if (index != -1) {
		if (index == 0) {
			if (val.length() == 1)
				throw new NumberFormatException("Zero length BigInteger");
			signum = -1;
			cursor = 1;
		} else {
			throw new NumberFormatException("Illegal embedded minus sign");
		}
	}
	
	// Skip leading zeros and compute number of digits in magnitude
	while (cursor < len &&
			Character.digit(val.charAt(cursor),radix) == 0)
		cursor++;
	if (cursor == len) {
		signum = 0;
		mag = ZERO.mag;
		return;
	} else {
		numDigits = len - cursor;
	}
	
	// Pre-allocate array of expected size. May be too large but can
	// never be too small. Typically exact.
	int numBits = (int)(((numDigits * bitsPerDigit[radix]) >>> 10) + 1);
	int numWords = (numBits + 31) /32;
	mag = new int[numWords];
	
	// Process first (potentially short) digit group
	int firstGroupLen = numDigits % digitsPerInt[radix];
	if (firstGroupLen == 0)
		firstGroupLen = digitsPerInt[radix];
	std::string group = val.substring(cursor, cursor += firstGroupLen);
	mag[mag.length - 1] = Integer.parseInt(group, radix);
	if (mag[mag.length - 1] < 0)
		throw new NumberFormatException("Illegal digit");
	
	// Process remaining digit groups
	int superRadix = intRadix[radix];
	int groupVal = 0;
	while (cursor < val.length()) {
		group = val.substring(cursor, cursor += digitsPerInt[radix]);
		groupVal = Integer.parseInt(group, radix);
		if (groupVal < 0)
			throw new NumberFormatException("Illegal digit");
		destructiveMulAdd(mag, superRadix, groupVal);
	}
	// Required for cases where the array was overallocated.
	mag = trustedStripLeadingZeroInts(mag);
}

BigInteger(char[] val):
	jcommon::Object()
{
	jcommon::Object::SetClassName("jmath::BigInteger");
	
	int cursor = 0, numDigits;
	int len = val.length;
	
	// Check for leading minus sign
	signum = 1;
	if (val[0] == '-') {
		if (len == 1)
			throw new NumberFormatException("Zero length BigInteger");
		signum = -1;
		cursor = 1;
	}
	
	// Skip leading zeros and compute number of digits in magnitude
	while (cursor < len && Character.digit(val[cursor], 10) == 0)
		cursor++;
	if (cursor == len) {
		signum = 0;
		mag = ZERO.mag;
		return;
	} else {
		numDigits = len - cursor;
	}
	
	// Pre-allocate array of expected size
	int numWords;
	if (len < 10) {
		numWords = 1;
	} else {    
		int numBits = (int)(((numDigits * bitsPerDigit[10]) >>> 10) + 1);
		numWords = (numBits + 31) /32;
	}
	mag = new int[numWords];
   	
	// Process first (potentially short) digit group
	int firstGroupLen = numDigits % digitsPerInt[10];
	if (firstGroupLen == 0)
		firstGroupLen = digitsPerInt[10];
	mag[mag.length-1] = parseInt(val, cursor,  cursor += firstGroupLen);
	
	// Process remaining digit groups
	while (cursor < len) {
		int groupVal = parseInt(val, cursor, cursor += digitsPerInt[10]);
		destructiveMulAdd(mag, intRadix[10], groupVal);
	}
	mag = trustedStripLeadingZeroInts(mag);
}

BigInteger::~BigInteger()
{
}

char * BigInteger::randomBits(int numBits, Random rnd) {
	if (numBits < 0)
		throw new IllegalArgumentException("numBits must be non-negative");
	int numBytes = (numBits+7)/8;
	char[] randomBits = new char[numBytes];
	
	// Generate random chars and mask out any excess bits
	if (numBytes > 0) {
		rnd.nextBytes(randomBits);
		int excessBits = 8*numBytes - numBits;
		randomBits[0] &= (1 << (8-excessBits)) - 1;
	}
	return randomBits;
}

BigInteger BigInteger::smallPrime(int bitLength, int certainty, Random rnd) {
	int magLen = (bitLength + 31) >>> 5;
	int temp[] = new int[magLen];
	int highBit = 1 << ((bitLength+31) & 0x1f);  // High bit of high int
	int highMask = (highBit << 1) - 1;  // Bits to keep in high int
	
	while(true) {
		// Construct a candidate
		for (int i=0; i<magLen; i++)
			temp[i] = rnd.nextInt();
		temp[0] = (temp[0] & highMask) | highBit;  // Ensure exact length
		if (bitLength > 2)
			temp[magLen-1] |= 1;  // Make odd if bitlen > 2
		
		BigInteger p = new BigInteger(temp, 1);
		
		// Do cheap "pre-test" if applicable
		if (bitLength > 6) {
			long r = p.remainder(SMALL_PRIME_PRODUCT).longValue();
			if ((r%3==0)  || (r%5==0)  || (r%7==0)  || (r%11==0) || 
					(r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) || 
					(r%29==0) || (r%31==0) || (r%37==0) || (r%41==0))
				continue; // Candidate is composite; try another
		}
		
		// All candidates of bitLength 2 and 3 are prime by this point
		if (bitLength < 4)
			return p;
		
		// Do expensive test if we survive pre-test (or it's inapplicable)
		if (p.primeToCertainty(certainty))
			return p;
	}
}

BigInteger BigInteger::largePrime(int bitLength, int certainty, Random rnd) {
	BigInteger p;
	p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
	p.mag[p.mag.length-1] &= 0xfffffffe;
	
	// Use a sieve length likely to contain the next prime number
	int searchLen = (bitLength / 20) * 64;
	BitSieve searchSieve = new BitSieve(p, searchLen);
	BigInteger candidate = searchSieve.retrieve(p, certainty);
	
	while ((candidate == NULL) || (candidate.bitLength() != bitLength)) {
		p = p.add(BigInteger.valueOf(2*searchLen));
		if (p.bitLength() != bitLength)
			p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
		p.mag[p.mag.length-1] &= 0xfffffffe;
		searchSieve = new BitSieve(p, searchLen);
		candidate = searchSieve.retrieve(p, certainty);
	}
	return candidate;
}

bool BigInteger::passesLucasLehmer() {
	BigInteger thisPlusOne = this.add(ONE);
	
	// Step 1
	int d = 5;
	while (jacobiSymbol(d, this) != -1) {
		// 5, -7, 9, -11, ...
		d = (d<0) ? Math.abs(d)+2 : -(d+2);
	}
	
	// Step 2
	BigInteger u = lucasLehmerSequence(d, thisPlusOne, this);
	
	// Step 3
	return u.mod(this).equals(ZERO);
}

int BigInteger::jacobiSymbol(int p, BigInteger n) {
	if (p == 0)
		return 0;
	
	// Algorithm and comments adapted from Colin Plumb's C library.
	int j = 1;
	int u = n.mag[n.mag.length-1];
	
	// Make p positive
	if (p < 0) {
		p = -p;
		int n8 = u & 7;
		if ((n8 == 3) || (n8 == 7))
			j = -j; // 3 (011) or 7 (111) mod 8
	}
	
	// Get rid of factors of 2 in p
	while ((p & 3) == 0)
		p >>= 2;
	if ((p & 1) == 0) {
		p >>= 1;
		if (((u ^ (u>>1)) & 2) != 0)
			j = -j;	// 3 (011) or 5 (101) mod 8
	}
	if (p == 1)
		return j;
	// Then, apply quadratic reciprocity
	if ((p & u & 2) != 0)	// p = u = 3 (mod 4)?
		j = -j;
	// And reduce u mod p
	u = n.mod(BigInteger.valueOf(p)).intValue();
	
	// Now compute Jacobi(u,p), u < p
	while (u != 0) {
		while ((u & 3) == 0)
			u >>= 2;
		if ((u & 1) == 0) {
			u >>= 1;
			if (((p ^ (p>>1)) & 2) != 0)
				j = -j;	// 3 (011) or 5 (101) mod 8
		}
		if (u == 1)
			return j;
		// Now both u and p are odd, so use quadratic reciprocity
		assert (u < p);
		int t = u; u = p; p = t;
		if ((u & p & 2) != 0) // u = p = 3 (mod 4)?
			j = -j;
		// Now u >= p, so it can be reduced
		u %= p;
	}
	return 0;
}

BigInteger BigInteger::lucasLehmerSequence(int z, BigInteger k, BigInteger n) {
	BigInteger d = BigInteger.valueOf(z);
	BigInteger u = ONE; BigInteger u2;
	BigInteger v = ONE; BigInteger v2;
	
	for (int i=k.bitLength()-2; i>=0; i--) {
		u2 = u.multiply(v).mod(n);
		
		v2 = v.square().add(d.multiply(u.square())).mod(n);
		if (v2.testBit(0)) {
			v2 = n.subtract(v2);
			v2.signum = - v2.signum;
		}
		v2 = v2.shiftRight(1);
		
		u = u2; v = v2;
		if (k.testBit(i)) {
			u2 = u.add(v).mod(n);
			if (u2.testBit(0)) {
				u2 = n.subtract(u2);
				u2.signum = - u2.signum;
			}
			u2 = u2.shiftRight(1);
			
			v2 = v.add(d.multiply(u)).mod(n);
			if (v2.testBit(0)) {
				v2 = n.subtract(v2);
				v2.signum = - v2.signum;
			}
			v2 = v2.shiftRight(1);
			
			u = u2; v = v2;
		}
	}
	return u;
}

bool BigInteger::passesMillerRabin(int iterations) {
	// Find a and m such that m is odd and this == 1 + 2**a * m
	BigInteger thisMinusOne = this.subtract(ONE);
	BigInteger m = thisMinusOne;
	int a = m.getLowestSetBit();
	m = m.shiftRight(a);
	
	// Do the tests
	Random rnd = new Random();
	for (int i=0; i<iterations; i++) {
		// Generate a uniform random on (1, this)
		BigInteger b;
		do {
			b = new BigInteger(this.bitLength(), rnd);
		} while (b.compareTo(ONE) <= 0 || b.compareTo(this) >= 0);
		
		int j = 0;
		BigInteger z = b.modPow(m, this);
		while(!((j==0 && z.equals(ONE)) || z.equals(thisMinusOne))) {
			if (j>0 && z.equals(ONE) || ++j==a)
				return false;
			z = z.modPow(TWO, this);
		}
	}
	return true;
}

BigInteger BigInteger::valueOf(int val[]) {
	return (val[0]>0 ? new BigInteger(val, 1) : new BigInteger(val));
}

int[] BigInteger::add(int[] x, int[] y) {
	// If x is shorter, swap the two arrays
	if (x.length < y.length) {
		int[] tmp = x;
		x = y;
		y = tmp;
	}
	
	int xIndex = x.length;
	int yIndex = y.length;
	int result[] = new int[xIndex];
	long sum = 0;
	
	// Add common parts of both numbers
	while(yIndex > 0) {
		sum = (x[--xIndex] & LONG_MASK) + 
			(y[--yIndex] & LONG_MASK) + (sum >>> 32);
		result[xIndex] = (int)sum;
	}
	
	// Copy remainder of longer number while carry propagation is required
	bool carry = (sum >>> 32 != 0);
	while (xIndex > 0 && carry)
		carry = ((result[--xIndex] = x[xIndex] + 1) == 0);
	
	// Copy remainder of longer number
	while (xIndex > 0)
		result[--xIndex] = x[xIndex];
	
	// Grow result if necessary
	if (carry) {
		int newLen = result.length + 1;
		int temp[] = new int[newLen];
		for (int i = 1; i<newLen; i++)
			temp[i] = result[i-1];
		temp[0] = 0x01;
		result = temp;
	}
	return result;
}

private int[] BigInteger::subtract(int[] big, int[] little) {
	int bigIndex = big.length;
	int result[] = new int[bigIndex];
	int littleIndex = little.length;
	long difference = 0;
	
	// Subtract common parts of both numbers
	while(littleIndex > 0) {
		difference = (big[--bigIndex] & LONG_MASK) - 
			(little[--littleIndex] & LONG_MASK) +
			(difference >> 32);
		result[bigIndex] = (int)difference;
	}
	
	// Subtract remainder of longer number while borrow propagates
	bool borrow = (difference >> 32 != 0);
	while (bigIndex > 0 && borrow)
		borrow = ((result[--bigIndex] = big[bigIndex] - 1) == -1);
	
	// Copy remainder of longer number
	while (bigIndex > 0)
		result[--bigIndex] = big[bigIndex];
	
	return result;
}

BigInteger BigInteger::multiply(BigInteger val) {
	if (signum == 0 || val.signum==0)
		return ZERO;
	
	int[] result = multiplyToLen(mag, mag.length, 
			val.mag, val.mag.length, NULL);
	result = trustedStripLeadingZeroInts(result);
	return new BigInteger(result, signum*val.signum);
}

int[] BigInteger::multiplyToLen(int[] x, int xlen, int[] y, int ylen, int[] z) {
	int xstart = xlen - 1;
	int ystart = ylen - 1;
	
	if (z == NULL || z.length < (xlen+ ylen))
		z = new int[xlen+ylen];
	
	long carry = 0;
	for (int j=ystart, k=ystart+1+xstart; j>=0; j--, k--) {
		long product = (y[j] & LONG_MASK) *
			(x[xstart] & LONG_MASK) + carry;
		z[k] = (int)product;
		carry = product >>> 32;
	}
	z[xstart] = (int)carry;
	
	for (int i = xstart-1; i >= 0; i--) {
		carry = 0;
		for (int j=ystart, k=ystart+1+i; j>=0; j--, k--) {
			long product = (y[j] & LONG_MASK) * 
				(x[i] & LONG_MASK) + 
				(z[k] & LONG_MASK) + carry;
			z[k] = (int)product;
			carry = product >>> 32;
		}
		z[i] = (int)carry;
	}
	return z;
}

	private BigInteger square() {
		if (signum == 0)
			return ZERO;
		int[] z = squareToLen(mag, mag.length, NULL);
		return new BigInteger(trustedStripLeadingZeroInts(z), 1);
	}

int[] BigInteger::squareToLen(int[] x, int len, int[] z) {
	int zlen = len << 1;
	if (z == NULL || z.length < zlen)
		z = new int[zlen];
	
	// Store the squares, right shifted one bit (i.e., divided by 2)
	int lastProductLowWord = 0;
	for (int j=0, i=0; j<len; j++) {
		long piece = (x[j] & LONG_MASK);
		long product = piece * piece;
		z[i++] = (lastProductLowWord << 31) | (int)(product >>> 33);
		z[i++] = (int)(product >>> 1);
		lastProductLowWord = (int)product;
	}
	
	// Add in off-diagonal sums
	for (int i=len, offset=1; i>0; i--, offset+=2) {
		int t = x[i-1];
		t = mulAdd(z, x, offset, i-1, t);
		addOne(z, offset-1, i, t);
	}
	
	// Shift back up and set low bit
	primitiveLeftShift(z, zlen, 1);
	z[zlen-1] |= x[len-1] & 1;
	
	return z;
}

int[] BigInteger::leftShift(int[] a, int len, int n) {
	int nInts = n >>> 5;
	int nBits = n&0x1F;
	int bitsInHighWord = bitLen(a[0]);
	
	// If shift can be done without recopy, do so
	if (n <= (32-bitsInHighWord)) {
		primitiveLeftShift(a, len, nBits);
		return a;
	} else { // Array must be resized
		if (nBits <= (32-bitsInHighWord)) {
			int result[] = new int[nInts+len];
			for (int i=0; i<len; i++)
				result[i] = a[i];
			primitiveLeftShift(result, result.length, nBits);
			return result;
		} else {
			int result[] = new int[nInts+len+1];
			for (int i=0; i<len; i++)
				result[i] = a[i];
			primitiveRightShift(result, result.length, 32 - nBits);
			return result;
		}
	}
}

void BigInteger::primitiveRightShift(int[] a, int len, int n) {
	int n2 = 32 - n;
	for (int i=len-1, c=a[i]; i>0; i--) {
		int b = c;
		c = a[i-1];
		a[i] = (c << n2) | (b >>> n);
	}
	a[0] >>>= n;
}

void BigInteger::primitiveLeftShift(int[] a, int len, int n) {
	if (len == 0 || n == 0)
		return;
	
	int n2 = 32 - n;
	for (int i=0, c=a[i], m=i+len-1; i<m; i++) {
		int b = c;
		c = a[i+1];
		a[i] = (b << n) | (c >>> n2);
	}
	a[len-1] <<= n;
}

int BigInteger::bitLength(int[] val, int len) {
	if (len==0)
		return 0;
	return ((len-1)<<5) + bitLen(val[0]);
}


int[] BigInteger::montReduce(int[] n, int[] mod, int mlen, int inv) {
	int c=0;
	int len = mlen;
	int offset=0;
	
	do {
		int nEnd = n[n.length-1-offset];
		int carry = mulAdd(n, mod, offset, mlen, inv * nEnd);
		c += addOne(n, offset, mlen, carry);
		offset++;
	} while(--len > 0);
	
	while(c>0)
		c += subN(n, mod, mlen);

        while (intArrayCmpToLen(n, mod, mlen) >= 0)
            subN(n, mod, mlen);

		return n;
}

int BigInteger::intArrayCmpToLen(int[] arg1, int[] arg2, int len) {
	for (int i=0; i<len; i++) {
		long b1 = arg1[i] & LONG_MASK;
		long b2 = arg2[i] & LONG_MASK;
		if (b1 < b2)
			return -1;
		if (b1 > b2)
			return 1;
	}
	return 0;
}

int BigInteger::subN(int[] a, int[] b, int len) {
	long sum = 0;
	
	while(--len >= 0) {
		sum = (a[len] & LONG_MASK) - 
			(b[len] & LONG_MASK) + (sum >> 32);
		a[len] = (int)sum;
	}
	
	return (int)(sum >> 32);
}

int BigInteger::mulAdd(int[] out, int[] in, int offset, int len, int k) {
	long kLong = k & LONG_MASK;
	long carry = 0;
	
	offset = out.length-offset - 1;
	for (int j=len-1; j >= 0; j--) {
		long product = (in[j] & LONG_MASK) * kLong +
			(out[offset] & LONG_MASK) + carry;
		out[offset--] = (int)product;
		carry = product >>> 32;
	}
	return (int)carry;
}

int BigInteger::addOne(int[] a, int offset, int mlen, int carry) {
	offset = a.length-1-mlen-offset;
	long t = (a[offset] & LONG_MASK) + (carry & LONG_MASK);
	
	a[offset] = (int)t;
	if ((t >>> 32) == 0)
		return 0;
	while (--mlen >= 0) {
		if (--offset < 0) { // Carry out of number
			return 1;
		} else {
			a[offset]++;
			if (a[offset] != 0)
				return 0;
		}
	}
	return 1;
}

BigInteger BigInteger::modPow2(BigInteger exponent, int p) {
	BigInteger result = valueOf(1);
	BigInteger baseToPow2 = this.mod2(p);
	int expOffset = 0;
	
	int limit = exponent.bitLength();
	
	if (this.testBit(0))
		limit = (p-1) < limit ? (p-1) : limit;
	
	while (expOffset < limit) {
		if (exponent.testBit(expOffset))
			result = result.multiply(baseToPow2).mod2(p);
		expOffset++;
		if (expOffset < limit)
			baseToPow2 = baseToPow2.square().mod2(p);
	}
	
	return result;
}

BigInteger 	BigInteger::mod2(int p) {
		if (bitLength() <= p)
			return this;
		
		// Copy remaining ints of mag
		int numInts = (p+31)/32;
		int[] mag = new int[numInts];
		for (int i=0; i<numInts; i++)
			mag[i] = this.mag[i + (this.mag.length - numInts)];
		
		// Mask out any excess bits
		int excessBits = (numInts << 5) - p;
		mag[0] &= (1L << (32-excessBits)) - 1;
		
		return (mag[0]==0 ? new BigInteger(1, mag) : new BigInteger(mag, 1));
	}

BigInteger BigInteger::oddModPow(BigInteger y, BigInteger z) {
	// Special case for exponent of one
	if (y.equals(ONE))
		return this;
	
	// Special case for base of zero
	if (signum==0)
		return ZERO;
	
	int[] base = (int[])mag.clone();
	int[] exp = y.mag;
	int[] mod = z.mag;
	int modLen = mod.length;
	
	// Select an appropriate window size
	int wbits = 0;
	int ebits = bitLength(exp, exp.length);
	// if exponent is 65537 (0x10001), use minimum window size
	if ((ebits != 17) || (exp[0] != 65537)) {
		while (ebits > bnExpModThreshTable[wbits]) {
			wbits++;
		}
	}
	
	// Calculate appropriate table size
	int tblmask = 1 << wbits;
	
	// Allocate table for precomputed odd powers of base in Montgomery form
	int[][] table = new int[tblmask][];
	for (int i=0; i<tblmask; i++)
		table[i] = new int[modLen];
	
	// Compute the modular inverse
	int inv = -MutableBigInteger.inverseMod32(mod[modLen-1]);
	
	// Convert base to Montgomery form
	int[] a = leftShift(base, base.length, modLen << 5);
	
	MutableBigInteger q = new MutableBigInteger(),
					  r = new MutableBigInteger(),
					  a2 = new MutableBigInteger(a),
					  b2 = new MutableBigInteger(mod);
	
	a2.divide(b2, q, r);
	table[0] = r.toIntArray();
	
	// Pad table[0] with leading zeros so its length is at least modLen
	if (table[0].length < modLen) {
		int offset = modLen - table[0].length;
		int[] t2 = new int[modLen];
		for (int i=0; i<table[0].length; i++)
			t2[i+offset] = table[0][i];
		table[0] = t2;
	}
	
	// Set b to the square of the base
	int[] b = squareToLen(table[0], modLen, NULL);
	b = montReduce(b, mod, modLen, inv);
	
	// Set t to high half of b
	int[] t = new int[modLen];
	for(int i=0; i<modLen; i++)
		t[i] = b[i];
	
	// Fill in the table with odd powers of the base        
	for (int i=1; i<tblmask; i++) {
		int[] prod = multiplyToLen(t, modLen, table[i-1], modLen, NULL);
		table[i] = montReduce(prod, mod, modLen, inv);
	}
	
	// Pre load the window that slides over the exponent
	int bitpos = 1 << ((ebits-1) & (32-1));
	
	int buf = 0;
	int elen = exp.length;
	int eIndex = 0;
	for (int i = 0; i <= wbits; i++) {
		buf = (buf << 1) | (((exp[eIndex] & bitpos) != 0)?1:0);
		bitpos >>>= 1;
		if (bitpos == 0) {
			eIndex++;
			bitpos = 1 << (32-1);
			elen--;
		}
	}
	
	int multpos = ebits;
	
	// The first iteration, which is hoisted out of the main loop
	ebits--;
	bool isone = true;
	
	multpos = ebits - wbits;
	while ((buf & 1) == 0) {
		buf >>>= 1;
		multpos++;
	}
	
	int[] mult = table[buf >>> 1];
	
	buf = 0;
	if (multpos == ebits)
		isone = false;
	
	// The main loop
	while(true) {
		ebits--;
		// Advance the window
		buf <<= 1;
		
		if (elen != 0) {
			buf |= ((exp[eIndex] & bitpos) != 0) ? 1 : 0;
			bitpos >>>= 1;
			if (bitpos == 0) {
				eIndex++;
				bitpos = 1 << (32-1);
				elen--;
			}
		}
		
		// Examine the window for pending multiplies
		if ((buf & tblmask) != 0) {
			multpos = ebits - wbits;
			while ((buf & 1) == 0) {
				buf >>>= 1;
				multpos++;
			}
			mult = table[buf >>> 1];
			buf = 0;
		}
		
		// Perform multiply
		if (ebits == multpos) {
			if (isone) {
				b = (int[])mult.clone();
				isone = false;
			} else {
				t = b;
				a = multiplyToLen(t, modLen, mult, modLen, a);
				a = montReduce(a, mod, modLen, inv);
				t = a; a = b; b = t;
			}
		}
		
		// Check if done
		if (ebits == 0)
			break;
		
		// Square the input
		if (!isone) {
			t = b;
			a = squareToLen(t, modLen, a);
			a = montReduce(a, mod, modLen, inv);
			t = a; a = b; b = t;
		}
	}

	// Convert result out of Montgomery form and return
	int[] t2 = new int[2*modLen];
	for(int i=0; i<modLen; i++)
		t2[i+modLen] = b[i];
	
	b = montReduce(t2, mod, modLen, inv);
	
	t2 = new int[modLen];
	for(int i=0; i<modLen; i++)
		t2[i] = b[i];
	
	return new BigInteger(1, t2);
}

int BigInteger::bitLen(int w) {
	// Binary search - decision tree (5 tests, rarely 6)
	return
		(w < 1<<15 ?
		 (w < 1<<7 ?
		  (w < 1<<3 ?
		   (w < 1<<1 ? (w < 1<<0 ? (w<0 ? 32 : 0) : 1) : (w < 1<<2 ? 2 : 3)) :
		   (w < 1<<5 ? (w < 1<<4 ? 4 : 5) : (w < 1<<6 ? 6 : 7))) :
		  (w < 1<<11 ?
		   (w < 1<<9 ? (w < 1<<8 ? 8 : 9) : (w < 1<<10 ? 10 : 11)) :
		   (w < 1<<13 ? (w < 1<<12 ? 12 : 13) : (w < 1<<14 ? 14 : 15)))) :
		 (w < 1<<23 ?
		  (w < 1<<19 ?
		   (w < 1<<17 ? (w < 1<<16 ? 16 : 17) : (w < 1<<18 ? 18 : 19)) :
		   (w < 1<<21 ? (w < 1<<20 ? 20 : 21) : (w < 1<<22 ? 22 : 23))) :
		  (w < 1<<27 ?
		   (w < 1<<25 ? (w < 1<<24 ? 24 : 25) : (w < 1<<26 ? 26 : 27)) :
		   (w < 1<<29 ? (w < 1<<28 ? 28 : 29) : (w < 1<<30 ? 30 : 31)))));
}

/* TODO::
// zero[i] is a string of i consecutive zeros. 
private std::string zeros[] = new std::string[64];
{
	zeros[63] =
		"000000000000000000000000000000000000000000000000000000000000000";
	for (int i=0; i<63; i++)
		zeros[i] = zeros[63].substring(0, i);
}
*/

int[] BigInteger::stripLeadingZeroInts(int val[]) {
	int charLength = val.length;
	int keep;
	
	// Find first nonzero char
	for (keep=0; keep<val.length && val[keep]==0; keep++)
		;
	
	int result[] = new int[val.length - keep];
	for(int i=0; i<val.length - keep; i++)
		result[i] = val[keep+i];
	
	return result;
}

int[] BigInteger::trustedStripLeadingZeroInts(int val[]) {
	int charLength = val.length;
	int keep;
	
	// Find first nonzero char
	for (keep=0; keep<val.length && val[keep]==0; keep++);
	
	// Only perform copy if necessary
	if (keep > 0) {
		int result[] = new int[val.length - keep];
		for(int i=0; i<val.length - keep; i++)
			result[i] = val[keep+i];
		return result; 
	}
	return val;
}

int[] BigInteger::stripLeadingZeroBytes(char a[]) {
	int charLength = a.length;
	int keep;
	
	// Find first nonzero char
	for (keep=0; keep<a.length && a[keep]==0; keep++)
		;
	
	// Allocate new array and copy relevant part of input array
	int intLength = ((charLength - keep) + 3)/4;
	int[] result = new int[intLength];
	int b = charLength - 1;
	for (int i = intLength-1; i >= 0; i--) {
		result[i] = a[b--] & 0xff;
		int charsRemaining = b - keep + 1;
		int charsToTransfer = Math.min(3, charsRemaining);
		for (int j=8; j <= 8*charsToTransfer; j += 8)
			result[i] |= ((a[b--] & 0xff) << j);
	}
	return result;
}

int[] BigInteger::makePositive(char a[]) {
	int keep, k;
	int charLength = a.length;
	
	// Find first non-sign (0xff) char of input
	for (keep=0; keep<charLength && a[keep]==-1; keep++)
		;
	
	
	/* Allocate output array.  If all non-sign chars are 0x00, we must
	 * allocate space for one extra output char. */
	for (k=keep; k<charLength && a[k]==0; k++)
		;
	
	int extraByte = (k==charLength) ? 1 : 0;
	int intLength = ((charLength - keep + extraByte) + 3)/4;
	int result[] = new int[intLength];
	
	/* Copy one's complement of input into output, leaving extra
	 * char (if it exists) == 0x00 */
	int b = charLength - 1;
	for (int i = intLength-1; i >= 0; i--) {
		result[i] = a[b--] & 0xff;
		int numBytesToTransfer = Math.min(3, b-keep+1);
		if (numBytesToTransfer < 0)
			numBytesToTransfer = 0;
		for (int j=8; j <= 8*numBytesToTransfer; j += 8)
			result[i] |= ((a[b--] & 0xff) << j);
		
		// Mask indicates which bits must be complemented
		int mask = -1 >>> (8*(3-numBytesToTransfer));
		result[i] = ~result[i] & mask;
	}
	
	// Add one to one's complement to generate two's complement
	for (int i=result.length-1; i>=0; i--) {
		result[i] = (int)((result[i] & LONG_MASK) + 1);
		if (result[i] != 0)
			break;
	}
	
	return result;
}

int[] BigInteger::makePositive(int a[]) {
	int keep, j;
	
	// Find first non-sign (0xffffffff) int of input
	for (keep=0; keep<a.length && a[keep]==-1; keep++)
		;
	
	/* Allocate output array.  If all non-sign ints are 0x00, we must
	 * allocate space for one extra output int. */
	for (j=keep; j<a.length && a[j]==0; j++)
		;
	int extraInt = (j==a.length ? 1 : 0);
	int result[] = new int[a.length - keep + extraInt];
	
	/* Copy one's complement of input into output, leaving extra
	 * int (if it exists) == 0x00 */
	for (int i = keep; i<a.length; i++)
		result[i - keep + extraInt] = ~a[i];
	
	// Add one to one's complement to generate two's complement
	for (int i=result.length-1; ++result[i]==0; i--)
		;
	
	return result;
}

int BigInteger::intLength() {
	return bitLength()/32 + 1;
}

/* Returns sign bit */
int BigInteger::signBit() {
	return (signum < 0 ? 1 : 0);
}

/* Returns an int of sign bits */
int BigInteger::signInt() {
	return (int) (signum < 0 ? -1 : 0);
}

int BigInteger::getInt(int n) {
	if (n < 0)
		return 0;
	if (n >= mag.length)
		return signInt();
	
	int magInt = mag[mag.length-n-1];
	
	return (int) (signum >= 0 ? magInt :
			(n <= firstNonzeroIntNum() ? -magInt : ~magInt));
}

int BigInteger::firstNonzeroIntNum() {
	if (firstNonzeroIntNum == -2) {
		// Search for the first nonzero int
		int i;
		for (i=mag.length-1; i>=0 && mag[i]==0; i--)
			;
		firstNonzeroIntNum = mag.length-i-1;
	}
	return firstNonzeroIntNum;
}

/**
 * Returns the mag array as an array of chars.
 */
char[] BigInteger::magSerializedForm() {
	int bitLen = (mag.length == 0 ? 0 :
			((mag.length - 1) << 5) + bitLen(mag[0]));
	int charLen = (bitLen + 7)/8;
	char[] result = new char[charLen];
	
	for (int i=charLen-1, charsCopied=4, intIndex=mag.length-1, nextInt=0;
			i>=0; i--) {
		if (charsCopied == 4) {
			nextInt = mag[intIndex--];
			charsCopied = 1;
		} else {
			nextInt >>>= 8;
			charsCopied++;
		}
		result[i] = (char)nextInt;
	}
	return result;
}

int BigInteger::parseInt(char[] source, int start, int end) {
	int result = Character.digit(source[start++], 10);
	if (result == -1)
		throw new NumberFormatException(new std::string(source));
	
	for (int index = start; index<end; index++) {
		int nextVal = Character.digit(source[index], 10);
		if (nextVal == -1)
			throw new NumberFormatException(new std::string(source));
		result = 10*result + nextVal;
	}
	
	return result;
}


void BigInteger::destructiveMulAdd(int[] x, int y, int z) {
	// Perform the multiplication word by word
	long ylong = y & LONG_MASK;
	long zlong = z & LONG_MASK;
	int len = x.length;
	
	long product = 0;
	long carry = 0;
	for (int i = len-1; i >= 0; i--) {
		product = ylong * (x[i] & LONG_MASK) + carry;
		x[i] = (int)product;
		carry = product >>> 32;
	}
	
	// Perform the addition
	long sum = (x[len-1] & LONG_MASK) + zlong;
	x[len-1] = (int)sum;
	carry = sum >>> 32;
	for (int i = len-2; i >= 0; i--) {
		sum = (x[i] & LONG_MASK) + carry;
		x[i] = (int)sum;
		carry = sum >>> 32;
	}
}

bool BigInteger::primeToCertainty(int certainty) {
	int rounds = 0;
	int n = (Math.min(certainty, Integer.MAX_VALUE-1)+1)/2;
	
	// The relationship between the certainty and the number of rounds
	// we perform is given in the draft standard ANSI X9.80, "PRIME
	// NUMBER GENERATION, PRIMALITY TESTING, AND PRIMALITY CERTIFICATES".
	int sizeInBits = this.bitLength();
	if (sizeInBits < 100) {
		rounds = 50;
		rounds = n < rounds ? n : rounds;
		return passesMillerRabin(rounds);
	}
	
	if (sizeInBits < 256) {
		rounds = 27;
	} else if (sizeInBits < 512) {
		rounds = 15;
	} else if (sizeInBits < 768) {
		rounds = 8;
	} else if (sizeInBits < 1024) {
		rounds = 4;
	} else {
		rounds = 2;
	}
	rounds = n < rounds ? n : rounds;
	
	return passesMillerRabin(rounds) && passesLucasLehmer();
}

int BigInteger::bitCnt(int val) {
	val -= (0xaaaaaaaa & val) >>> 1;
	val = (val & 0x33333333) + ((val >>> 2) & 0x33333333);
	val = val + (val >>> 4) & 0x0f0f0f0f;
	val += val >>> 8;
	val += val >>> 16;
	return val & 0xff;
}

int BigInteger::trailingZeroCnt(int val) {
	// Loop unrolled for performance
	int charVal = val & 0xff;
	if (charVal != 0)
		return trailingZeroTable[charVal];
	
	charVal = (val >>> 8) & 0xff;
	if (charVal != 0)
		return trailingZeroTable[charVal] + 8;
	
	charVal = (val >>> 16) & 0xff;
	if (charVal != 0)
		return trailingZeroTable[charVal] + 16;
	
	charVal = (val >>> 24) & 0xff;
	return trailingZeroTable[charVal] + 24;
}


int BigInteger::bitCnt(int val) {
	val -= (0xaaaaaaaa & val) >>> 1;
	val = (val & 0x33333333) + ((val >>> 2) & 0x33333333);
	val = val + (val >>> 4) & 0x0f0f0f0f;
	val += val >>> 8;
	val += val >>> 16;
	return val & 0xff;
}

int BigInteger::trailingZeroCnt(int val) {
	// Loop unrolled for performance
	int charVal = val & 0xff;
	if (charVal != 0)
		return trailingZeroTable[charVal];
	
	charVal = (val >>> 8) & 0xff;
	if (charVal != 0)
		return trailingZeroTable[charVal] + 8;
	
	charVal = (val >>> 16) & 0xff;
	if (charVal != 0)
		return trailingZeroTable[charVal] + 16;
	
	charVal = (val >>> 24) & 0xff;
	return trailingZeroTable[charVal] + 24;
}

int BigInteger::intArrayCmp(int[] arg1, int[] arg2) {
	if (arg1.length < arg2.length)
		return -1;
	if (arg1.length > arg2.length)
		return 1;
	
	// Argument lengths are equal; compare the values
	for (int i=0; i<arg1.length; i++) {
		long b1 = arg1[i] & LONG_MASK;
		long b2 = arg2[i] & LONG_MASK;
		if (b1 < b2)
			return -1;
		if (b1 > b2)
			return 1;
	}
	return 0;
}

BigInteger BigInteger::probablePrime(int bitLength, Random rnd) {
	if (bitLength < 2)
		throw new ArithmeticException("bitLength < 2");
	
	// The cutoff of 95 was chosen empirically for best performance
	return (bitLength < SMALL_PRIME_THRESHOLD ?
			smallPrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd) :
			largePrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd));
}

BigInteger BigInteger::nextProbablePrime() {
	if (this.signum < 0)
		throw new ArithmeticException("start < 0: " + this);
	
	// Handle trivial cases
	if ((this.signum == 0) || this.equals(ONE))
		return TWO;
	
	BigInteger result = this.add(ONE);
	
	// Fastpath for small numbers
	if (result.bitLength() < SMALL_PRIME_THRESHOLD) {
	   	
   		// Ensure an odd number
		if (!result.testBit(0))
			result = result.add(ONE);
		
		while(true) {
			// Do cheap "pre-test" if applicable
			if (result.bitLength() > 6) {
				long r = result.remainder(SMALL_PRIME_PRODUCT).longValue();
				if ((r%3==0)  || (r%5==0)  || (r%7==0)  || (r%11==0) || 
						(r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) || 
						(r%29==0) || (r%31==0) || (r%37==0) || (r%41==0)) {
					result = result.add(TWO);
					continue; // Candidate is composite; try another
				}
			}
			
			// All candidates of bitLength 2 and 3 are prime by this point
			if (result.bitLength() < 4)
				return result;
			
			// The expensive test
			if (result.primeToCertainty(DEFAULT_PRIME_CERTAINTY))
				return result;
			
			result = result.add(TWO);
		}
	}
	
	// Start at previous even number
	if (result.testBit(0))
		result = result.subtract(ONE);
	
	// Looking for the next large prime
	int searchLen = (result.bitLength() / 20) * 64;
	
	while(true) {
		BitSieve searchSieve = new BitSieve(result, searchLen);
		BigInteger candidate = searchSieve.retrieve(result,
				DEFAULT_PRIME_CERTAINTY);
		if (candidate != NULL)
			return candidate;
		result = result.add(BigInteger.valueOf(2 * searchLen));
	}
}

BigInteger BigInteger::valueOf(long val) {
	// If -MAX_CONSTANT < val < MAX_CONSTANT, return stashed constant
	if (val == 0)
		return ZERO;
	if (val > 0 && val <= MAX_CONSTANT)
		return posConst[(int) val];
	else if (val < 0 && val >= -MAX_CONSTANT)
		return negConst[(int) -val];
	
	return new BigInteger(val);
}

BigInteger BigInteger::add(BigInteger val) {
	int[] resultMag;
	if (val.signum == 0)
		return this;
	if (signum == 0)
		return val;
	if (val.signum == signum)
		return new BigInteger(add(mag, val.mag), signum);
	
	int cmp = intArrayCmp(mag, val.mag);
	if (cmp==0)
		return ZERO;
	resultMag = (cmp>0 ? subtract(mag, val.mag)
			: subtract(val.mag, mag));
	resultMag = trustedStripLeadingZeroInts(resultMag);
	
	return new BigInteger(resultMag, cmp*signum);
}

BigInteger BigInteger::divide(BigInteger val) {
	MutableBigInteger q = new MutableBigInteger(),
					  r = new MutableBigInteger(),
					  a = new MutableBigInteger(this.mag),
					  b = new MutableBigInteger(val.mag);
	
	a.divide(b, q, r);
	return new BigInteger(q, this.signum * val.signum);
}

BigInteger[] BigInteger::divideAndRemainder(BigInteger val) {
	BigInteger[] result = new BigInteger[2];
	MutableBigInteger q = new MutableBigInteger(),
					  r = new MutableBigInteger(),
					  a = new MutableBigInteger(this.mag),
					  b = new MutableBigInteger(val.mag);
	a.divide(b, q, r);
	result[0] = new BigInteger(q, this.signum * val.signum);
	result[1] = new BigInteger(r, this.signum);
	return result;
}

BigInteger BigInteger::remainder(BigInteger val) {
	MutableBigInteger q = new MutableBigInteger(),
					  r = new MutableBigInteger(),
					  a = new MutableBigInteger(this.mag),
					  b = new MutableBigInteger(val.mag);
	
	a.divide(b, q, r);
	return new BigInteger(r, this.signum);
}

BigInteger BigInteger::pow(int exponent) {
	if (exponent < 0)
		throw new ArithmeticException("Negative exponent");
	if (signum==0)
		return (exponent==0 ? ONE : this);
	
	// Perform exponentiation using repeated squaring trick
	int newSign = (signum<0 && (exponent&1)==1 ? -1 : 1);
	int[] baseToPow2 = this.mag;
	int[] result = {1};
	
	while (exponent != 0) {
		if ((exponent & 1)==1) {
			result = multiplyToLen(result, result.length, 
					baseToPow2, baseToPow2.length, NULL);
			result = trustedStripLeadingZeroInts(result);
		}
		if ((exponent >>>= 1) != 0) {
			baseToPow2 = squareToLen(baseToPow2, baseToPow2.length, NULL);
			baseToPow2 = trustedStripLeadingZeroInts(baseToPow2);
		}
	}
	return new BigInteger(result, newSign);
}

BigInteger BigInteger::gcd(BigInteger val) {
	if (val.signum == 0)
		return this.abs();
	else if (this.signum == 0)
		return val.abs();
	
	MutableBigInteger a = new MutableBigInteger(this);
	MutableBigInteger b = new MutableBigInteger(val);
	
	MutableBigInteger result = a.hybridGCD(b);
	
	return new BigInteger(result, 1);
}

BigInteger BigInteger::abs() {
	return (signum >= 0 ? this : this.negate());
}

BigInteger BigInteger::negate() {
	return new BigInteger(this.mag, -this.signum);
}

int BigInteger::signum() {
	return this.signum;
}

BigInteger BigInteger::mod(BigInteger m) {
	if (m.signum <= 0)
		throw new ArithmeticException("BigInteger: modulus not positive");
	
	BigInteger result = this.remainder(m);
	return (result.signum >= 0 ? result : result.add(m));
}

BigInteger BigInteger::modPow(BigInteger exponent, BigInteger m) {
	if (m.signum <= 0)
		throw new ArithmeticException("BigInteger: modulus not positive");
	
	// Trivial cases
	if (exponent.signum == 0)
		return (m.equals(ONE) ? ZERO : ONE);
	
	if (this.equals(ONE))
		return (m.equals(ONE) ? ZERO : ONE);
	
	if (this.equals(ZERO) && exponent.signum >= 0)
		return ZERO;
	
	if (this.equals(negConst[1]) && (!exponent.testBit(0)))
		return (m.equals(ONE) ? ZERO : ONE);
	
	bool invertResult;
	if ((invertResult = (exponent.signum < 0)))
		exponent = exponent.negate();
	
	BigInteger base = (this.signum < 0 || this.compareTo(m) >= 0
 			? this.mod(m) : this);
	BigInteger result;
	if (m.testBit(0)) { // odd modulus
		result = base.oddModPow(exponent, m);
	} else {
		/*
		 * Even modulus.  Tear it into an "odd part" (m1) and power of two
		 * (m2), exponentiate mod m1, manually exponentiate mod m2, and
		 * use Chinese Remainder Theorem to combine results.
		 */
		
		// Tear m apart into odd part (m1) and power of 2 (m2)
		int p = m.getLowestSetBit();   // Max pow of 2 that divides m
		
		BigInteger m1 = m.shiftRight(p);  // m/2**p
		BigInteger m2 = ONE.shiftLeft(p); // 2**p
		
		// Calculate new base from m1
		BigInteger base2 = (this.signum < 0 || this.compareTo(m1) >= 0
				? this.mod(m1) : this);
		
		// Caculate (base ** exponent) mod m1.
		BigInteger a1 = (m1.equals(ONE) ? ZERO :
				base2.oddModPow(exponent, m1));
		
		// Calculate (this ** exponent) mod m2
		BigInteger a2 = base.modPow2(exponent, p);
		
		// Combine results using Chinese Remainder Theorem
		BigInteger y1 = m2.modInverse(m1);
		BigInteger y2 = m1.modInverse(m2);
		
		result = a1.multiply(m2).multiply(y1).add
	   		(a2.multiply(m1).multiply(y2)).mod(m);
	}
	
	return (invertResult ? result.modInverse(m) : result);
}

BigInteger BigInteger::modInverse(BigInteger m) {
	if (m.signum != 1)
		throw new ArithmeticException("BigInteger: modulus not positive");
	
	if (m.equals(ONE))
		return ZERO;
	
	// Calculate (this mod m)
	BigInteger modVal = this;
	if (signum < 0 || (intArrayCmp(mag, m.mag) >= 0))
		modVal = this.mod(m);
	
	if (modVal.equals(ONE))
		return ONE;
	
	MutableBigInteger a = new MutableBigInteger(modVal);
	MutableBigInteger b = new MutableBigInteger(m);
  	
	MutableBigInteger result = a.mutableModInverse(b);  
	return new BigInteger(result, 1);
}

BigInteger BigInteger::shiftLeft(int n) {
	if (signum == 0)
		return ZERO;
	if (n==0)
		return this;
	if (n<0)
		return shiftRight(-n);
	
	int nInts = n >>> 5;
	int nBits = n & 0x1f;
	int magLen = mag.length;
	int newMag[] = NULL;
	
	if (nBits == 0) {
		newMag = new int[magLen + nInts];
		for (int i=0; i<magLen; i++)
			newMag[i] = mag[i];
	} else {
		int i = 0;
		int nBits2 = 32 - nBits;
		int highBits = mag[0] >>> nBits2;
		if (highBits != 0) {
			newMag = new int[magLen + nInts + 1];
			newMag[i++] = highBits;
		} else {
			newMag = new int[magLen + nInts];
		}
		int j=0;
		while (j < magLen-1)
			newMag[i++] = mag[j++] << nBits | mag[j] >>> nBits2;
		newMag[i] = mag[j] << nBits;
	}
	
	return new BigInteger(newMag, signum);
}

BigInteger BigInteger::shiftRight(int n) {
	if (n==0)
		return this;
	if (n<0)
		return shiftLeft(-n);
	
	int nInts = n >>> 5;
	int nBits = n & 0x1f;
	int magLen = mag.length;
	int newMag[] = NULL;
	
	// Special case: entire contents shifted off the end
	if (nInts >= magLen)
		return (signum >= 0 ? ZERO : negConst[1]);
	
	if (nBits == 0) {
		int newMagLen = magLen - nInts;
		newMag = new int[newMagLen];
		for (int i=0; i<newMagLen; i++)
			newMag[i] = mag[i];
	} else {
		int i = 0;
		int highBits = mag[0] >>> nBits;
		if (highBits != 0) {
			newMag = new int[magLen - nInts];
			newMag[i++] = highBits;
		} else {
			newMag = new int[magLen - nInts -1];
		}
		
		int nBits2 = 32 - nBits;
		int j=0;
		while (j < magLen - nInts - 1)
			newMag[i++] = (mag[j++] << nBits2) | (mag[j] >>> nBits);
	}
	
	if (signum < 0) {
		// Find out whether any one-bits were shifted off the end.
		bool onesLost = false;
		for (int i=magLen-1, j=magLen-nInts; i>=j && !onesLost; i--)
			onesLost = (mag[i] != 0);
		if (!onesLost && nBits != 0)
			onesLost = (mag[magLen - nInts - 1] << (32 - nBits) != 0);
		
		if (onesLost)
			newMag = javaIncrement(newMag);
        }
	
	return new BigInteger(newMag, signum);
}

int[] BigInteger::javaIncrement(int[] val) {
	bool done = false;
	int lastSum = 0;
	for (int i=val.length-1;  i >= 0 && lastSum == 0; i--)
		lastSum = (val[i] += 1);
	if (lastSum == 0) {
		val = new int[val.length+1];
		val[0] = 1;
	}
	return val;
}

BigInteger BigInteger::and(BigInteger val) {
	int[] result = new int[Math.max(intLength(), val.intLength())];
	for (int i=0; i<result.length; i++)
		result[i] = (int) (getInt(result.length-i-1)
				& val.getInt(result.length-i-1));
	
	return valueOf(result);
}

BigInteger BigInteger::or(BigInteger val) {
	int[] result = new int[Math.max(intLength(), val.intLength())];
	for (int i=0; i<result.length; i++)
		result[i] = (int) (getInt(result.length-i-1)
				| val.getInt(result.length-i-1));
	
	return valueOf(result);
}

BigInteger BigInteger::xor(BigInteger val) {
	int[] result = new int[Math.max(intLength(), val.intLength())];
	for (int i=0; i<result.length; i++)
		result[i] = (int) (getInt(result.length-i-1)
				^ val.getInt(result.length-i-1));
	
	return valueOf(result);
}

BigInteger BigInteger::not() {
	int[] result = new int[intLength()];
	for (int i=0; i<result.length; i++)
		result[i] = (int) ~getInt(result.length-i-1);
	
	return valueOf(result);
}

BigInteger BigInteger::andNot(BigInteger val) {
	int[] result = new int[Math.max(intLength(), val.intLength())];
	for (int i=0; i<result.length; i++)
		result[i] = (int) (getInt(result.length-i-1)
				& ~val.getInt(result.length-i-1));
	
	return valueOf(result);
}

bool BigInteger::testBit(int n) {
	if (n<0)
		throw new ArithmeticException("Negative bit address");
	
	return (getInt(n/32) & (1 << (n%32))) != 0;
}

BigInteger BigInteger::setBit(int n) {
	if (n<0)
		throw new ArithmeticException("Negative bit address");
	
	int intNum = n/32;
	int[] result = new int[Math.max(intLength(), intNum+2)];
	
	for (int i=0; i<result.length; i++)
		result[result.length-i-1] = getInt(i);
	
	result[result.length-intNum-1] |= (1 << (n%32));
	
	return valueOf(result);
}

BigInteger BigInteger::clearBit(int n) {
	if (n<0)
		throw new ArithmeticException("Negative bit address");
	
	int intNum = n/32;
	int[] result = new int[Math.max(intLength(), (n+1)/32+1)];
	
	for (int i=0; i<result.length; i++)
		result[result.length-i-1] = getInt(i);
	
	result[result.length-intNum-1] &= ~(1 << (n%32));
	
	return valueOf(result);
}

BigInteger BigInteger::flipBit(int n) {
	if (n<0)
		throw new ArithmeticException("Negative bit address");
	
	int intNum = n/32;
	int[] result = new int[Math.max(intLength(), intNum+2)];
	
	for (int i=0; i<result.length; i++)
		result[result.length-i-1] = getInt(i);
	
	result[result.length-intNum-1] ^= (1 << (n%32));
	
	return valueOf(result);
}

int BigInteger::getLowestSetBit() {
	/*
	 * Initialize lowestSetBit field the first time this method is
	 * executed. This method depends on the atomicity of int modifies;
	 * without this guarantee, it would have to be synchronized.
	 */
	if (lowestSetBit == -2) {
		if (signum == 0) {
			lowestSetBit = -1;
		} else {
			// Search for lowest order nonzero int
			int i,b;
			for (i=0; (b = getInt(i))==0; i++)
				;
			lowestSetBit = (i << 5) + trailingZeroCnt(b);
		}
	}
	return lowestSetBit;
}

int BigInteger::bitLength() {
	/*
	 * \brief Initialize bitLength field the first time this method is executed.
	 * This method depends on the atomicity of int modifies; without
	 * this guarantee, it would have to be synchronized.
	 */
	if (bitLength == -1) {
		if (signum == 0) {
			bitLength = 0;
		} else {
			// Calculate the bit length of the magnitude
			int magBitLength = ((mag.length-1) << 5) + bitLen(mag[0]);
			
			if (signum < 0) {
				// Check if magnitude is a power of two
				bool pow2 = (bitCnt(mag[0]) == 1);
				for(int i=1; i<mag.length && pow2; i++)
					pow2 = (mag[i]==0);
				
				bitLength = (pow2 ? magBitLength-1 : magBitLength);
			} else {
				bitLength = magBitLength;
			}
		}
	}
	return bitLength;
}

int BigInteger::bitCount() {
	/*
	 * Initialize bitCount field the first time this method is executed.
	 * This method depends on the atomicity of int modifies; without
	 * this guarantee, it would have to be synchronized.
	 */
	if (bitCount == -1) {
		// Count the bits in the magnitude
		int magBitCount = 0;
		for (int i=0; i<mag.length; i++)
			magBitCount += bitCnt(mag[i]);
		
		if (signum < 0) {
			// Count the trailing zeros in the magnitude
			int magTrailingZeroCount = 0, j;
			for (j=mag.length-1; mag[j]==0; j--)
				magTrailingZeroCount += 32;
			magTrailingZeroCount +=
				trailingZeroCnt(mag[j]);
			
			bitCount = magBitCount + magTrailingZeroCount - 1;
		} else {
			bitCount = magBitCount;
		}
	}
	return bitCount;
}

BigInteger BigInteger::subtract(BigInteger val) {
	int[] resultMag;
	if (val.signum == 0)
		return this;
	if (signum == 0)
		return val.negate();
	if (val.signum != signum)
		return new BigInteger(add(mag, val.mag), signum);
	
	int cmp = intArrayCmp(mag, val.mag);
	if (cmp==0)
		return ZERO;
	resultMag = (cmp>0 ? subtract(mag, val.mag)
			: subtract(val.mag, mag));
	resultMag = trustedStripLeadingZeroInts(resultMag);
	return new BigInteger(resultMag, cmp*signum);
}

bool BigInteger::isProbablePrime(int certainty) {
	if (certainty <= 0)
		return true;
	BigInteger w = this.abs();
	if (w.equals(TWO))
		return true;
	if (!w.testBit(0) || w.equals(ONE))
		return false;
	
	return w.primeToCertainty(certainty);
}

int BigInteger::compareTo(BigInteger val) {
	return (signum==val.signum
			? signum*intArrayCmp(mag, val.mag)
			: (signum>val.signum ? 1 : -1));
}

bool BigInteger::equals(Object x) {
	// This test is just an optimization, which may or may not help
	if (x == this)
		return true;
	
	if (!(x instanceof BigInteger))
		return false;
	BigInteger xInt = (BigInteger) x;
	
	if (xInt.signum != signum || xInt.mag.length != mag.length)
		return false;
	
	for (int i=0; i<mag.length; i++)
		if (xInt.mag[i] != mag[i])
			return false;
	
	return true;
}

BigInteger BigInteger::min(BigInteger val) {
	return (compareTo(val)<0 ? this : val);
}

BigInteger BigInteger::max(BigInteger val) {
	return (compareTo(val)>0 ? this : val);
}

int BigInteger::hashCode() {
	int hashCode = 0;
	
	for (int i=0; i<mag.length; i++)
		hashCode = (int)(31*hashCode + (mag[i] & LONG_MASK));
	
	return hashCode * signum;
}

std::string BigInteger::string(int radix) {
	if (signum == 0)
		return "0";
	if (radix < Character.MIN_RADIX || radix > Character.MAX_RADIX)
		radix = 10;
	
	// Compute upper bound on number of digit groups and allocate space
	int maxNumDigitGroups = (4*mag.length + 6)/7;
	std::string digitGroup[] = new std::string[maxNumDigitGroups];
	
	// Translate number to string, a digit group at a time
	BigInteger tmp = this.abs();
	int numGroups = 0;
	while (tmp.signum != 0) {
		BigInteger d = longRadix[radix];
		
		MutableBigInteger q = new MutableBigInteger(),
						  r = new MutableBigInteger(),
						  a = new MutableBigInteger(tmp.mag),
						  b = new MutableBigInteger(d.mag);
		a.divide(b, q, r);
		BigInteger q2 = new BigInteger(q, tmp.signum * d.signum);
		BigInteger r2 = new BigInteger(r, tmp.signum * d.signum);
		
		digitGroup[numGroups++] = Long.tostd::string(r2.longValue(), radix);
		tmp = q2;
	}
	
	// Put sign (if any) and first digit group into result buffer
	std::stringBuilder buf = new std::stringBuilder(numGroups*digitsPerLong[radix]+1);
	if (signum<0)
		buf.append('-');
	buf.append(digitGroup[numGroups-1]);
	
	// Append remaining digit groups padded with leading zeros
	for (int i=numGroups-2; i>=0; i--) {
		// Prepend (any) leading zeros for this digit group
		int numLeadingZeros = digitsPerLong[radix]-digitGroup[i].length();
		if (numLeadingZeros != 0)
			buf.append(zeros[numLeadingZeros]);
		buf.append(digitGroup[i]);
	}
	return buf.tostd::string();
}

std::string BigInteger::string() {
	return tostd::string(10);
}

char[] BigInteger::toByteArray() {
	int charLen = bitLength()/8 + 1;
	char[] charArray = new char[charLen];
	
	for (int i=charLen-1, charsCopied=4, nextInt=0, intIndex=0; i>=0; i--) {
		if (charsCopied == 4) {
			nextInt = getInt(intIndex++);
			charsCopied = 1;
		} else {
			nextInt >>>= 8;
			charsCopied++;
		}
		charArray[i] = (char)nextInt;
	}
	return charArray;
}

int BigInteger::intValue() {
	int result = 0;
	result = getInt(0);
	return result;
}

long BigInteger::longValue() {
	long result = 0;
	
	for (int i=1; i>=0; i--)
		result = (result << 32) + (getInt(i) & LONG_MASK);
	return result;
}

float BigInteger::floatValue() {
	// Somewhat inefficient, but guaranteed to work.
	return Float.valueOf(this.tostd::string()).floatValue();
}

double BigInteger::doubleValue() {
	// Somewhat inefficient, but guaranteed to work.
	return Double.valueOf(this.tostd::string()).doubleValue();
}

};