// Emacs settings: -*- mode: Fundamental; tab-width: 4; -*-

////////////////////////////////////////////////////////////////////////////
//                                                                        //
// Javascript arbitrary precision integers: Bignum                        //
//                                                                        //
// Copyright (c) 2009-2010, Andrew Birrell                                //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

// This library implements arbitrary-precision integers in Javascript,
// including support for negative integers and NaN.  It includes basic
// arithmetic, and a few number theory operations (GCD, modular inverse,
// and primality testing).
//
// The public entry points are as follows.  See the source for more detail.
//
//    new Bignum(x) ... create BigNum object from a Number or decimal String
//    bigZero       ... = new Bignum(0)
//    bigOne        ... = new Bignum(1)
//    bigTwo        ... = new Bignum(2)
//    a.toString()  ... a decimal string for "a"
//    a.absComp(b)  ... (-1, 0, +1) if |a| is (lt, eq, gt) |b|
//    a.comp(b)     ... (-1, 0, +1) if a is (lt, eq, gt) b
//    a.negate()    ... -a
//    a.add(b)      ... a + b
//    a.sub(b)      ... a - b
//    a.mul(b)      ... a * b
//    a.mul(b, m)   ... (a * b) MOD m
//    a.mulInt(n)   ... a * n for ordinary positive integer n
//    a.div(b)      ... a / b, rounding towards 0
//    a.mod(b)      ... a MOD b, s.t. a = a * (a / b) + a MOD b
//    a.pow(b)      ... a ^ b
//    a.pow(b, m)   ... (a ^ b) MOD m
//    a.gcd(b)      ... gcd(|a|, |b|)
//    a.inverse(b)  ... c s.t. (a * c) MOD b = 1; requires gcd(a, b) = 1
//    a.isPrime()   ... true iff a is probably prime (to 1 in 2^80)
//    a.randomLT()  ... random number in [0 .. |a|-1]
//    bigRandom(n)  ... random number in [0 .. 2^n-1]
//    bigPrime(n)   ... random probable prime in [2^(n-1) .. 2^n-1]
//
// In general, our semantics mirror the fine print of EcmaScript Number
// semantics, except for our omission of infinity and negative zero:
//
//   - errors are convert into "NaN"; the "NaN" flag is a relevant string;
//   - "NaN" propagates through most operations ("pow(0)" is a bit special);
//   - the Bignum constructor behaves like "Number(x)" for all x;
//   - division rounds towards zero;
//   - remainder = dividend - quotient * divisor.
//
// "pow" limits the allowed results to terminate infeasible computation.
//
// Bignums are represented as objects containing a vector of positive
// integers, each integer < radix, together with a sign bit and a "NaN"
// flag.  The objects are immutable, except while being created during an
// operation.
//
// Overall, the performance is moderately good, for Javascript.  Multiply
// uses the Karatsuba algorithm for numbers large enough for it to be
// faster.  Divmod uses classic schoolroom long division.  For example,
// "a.pow(511)" for 512-bit "a" takes about 400 msec on a 2.4 GHz Core Duo
// Mac with Safari 5.
//
// The FFT multiplication algorithm is asymptotically faster than Karatsuba,
// but it wouldn't help until very much larger numbers.  Neither Newton-
// Raphson division nor Montgomery reduction (for a^b mod m) help unless
// you have much faster multiply.  The primary candidate for optimizing the
// division operations with the existing algorithms is the inner loop of
// bigMulIntSub.
//
// This library is provided, if at all, entirely "as is"; use it at your
// own risk.  Don't even dream of using the random numbers or the prime
// generator to protect anything of value.
//
// The library is copyright: if you want to use it, ask first.


var radix = 64*1024*1024;
	// As coded below, we require exact representation of integers up to
	// 2*(radix*radix-1) ... see the inner loop of "bigMulClassic".
	//
	// EcmaScript represents numbers as 64-bit IEEE floats, with exact
	// integers up to 2^53 (roughly 9 * 10^15).  That allows radix 2^26.
	//
	// (With Safari 4.0.4 on MacOS 10.5 on PowerPC radix > 2^20 fails).
	//
	// "bigPowMS" assumes radix is a power of 2.
	//
var radixBits = 26;
	// radix <= 2 ^ radixBits; for Miller-Rabin.

function Bignum(n) {
	// Constructor for a Bignum with given value.
	//
	this.neg = false;
	if (n === 0 || n === false || n === null || n === "") {
		// Values that map to 0, done first for efficiency
		var v = new Array();
		v.push(0);
		this.v = v;
	} else if (!n) {
		// Values that map to NaN: the number "NaN" and undefined
		this.NaN = true;
	} else if (typeof(n) == "number" || typeof(n) == "boolean") {
		// Native number, or boolean "true" converted implicitly.
		if (n < 0) {
			this.neg = true;
			n = -n;
		}
		var v = new Array();
		while (n > 0) {
			var n2 = Math.floor(n / radix);
			v.push(n - n2 * radix);
			n = n2;
		}
		if (v.length == 0) v.push(0);
		this.v = v;
	} else {
		// String, Object, Function, and anything not in EcmaScript spec.
		var temp = null;
		if (typeof(n) == "string") {
			while (n.length > 0 && n.charAt(0) == " ") {
				n = n.substr(1,n.length);
			}
			if (n.length > 0 && n.charAt(0) == "-") {
				n = n.substr(1, n.length);
				this.neg = true;
			}
			for (var i = 0; i < n.length; i++) {
				var c = n.charAt(i);
				if (c != " ") {
					var next = new Bignum(Number(c));
					if (next.NaN) {
						temp = null;
						break;
					}
					temp = (temp ? temp.mulInt(10).add(next) : next);
				}
			}
			if (!temp) {
				this.NaN = n;
			} else {
				this.v = temp.v;
			}
		} else {
			this.NaN = "syntax error in constant";
		}
	}
}

var bigZero = new Bignum(0);
var bigOne = new Bignum(1);
var bigTwo = new Bignum(2);

function bigClone(a) {
	// INTERNAL: returns a copy of a's integer array.
	//
	return a.v.slice(0);
}

function bigDivInt(v, n) {
	// INTERNAL: used to convert a Bignum to binary (for "pow") or decimal
	// (for "toString")
	//
	// Divide the Bignum that would be represented by "v" by the regular
	// positive integer "n", destructively, and return the remainder.
	//
	// Does not remove leading zeros in "v"; the caller should.
	//
	var limit = v.length-1;
	var rem = 0;
	for (var i = limit; i >= 0; i--) {
		var x = rem * radix + v[i];
		var quotient = Math.floor(x / n);
		rem = x - quotient * n;
		v[i] = quotient;
	}
	return rem;
}

Bignum.prototype.toString = function() {
	// Convert "this" to a string of decimal digits, or NaN.
	//
	if (this.NaN) {
		return "NaN" + (typeof(this.NaN) == "string" ?
			" (" + this.NaN + ")" : "");
	}
	var clump = 10000000;
	var digits = 7;
	var t = bigClone(this);
	var res = new Array();
	while (t.length > 0) {
		var limit = t.length-1;
		if (t[limit] == 0) {
			t.length = limit;
		} else {
			res.push(bigDivInt(t, clump));
		}
	}
	if (res.length == 0) return "0";
	var s = "";
	for (var i = 0; i < res.length; i++) {
		// "res" has l.s. digits first
		var d = res[i].toString();
		if (i < res.length-1) {
			// insert leading zeros, except at  the very front
			var dLen = d.length;
			for (var j = 0; j < digits - dLen; j++) {
				d = "0" + d;
			}
		}
		s = d + s;
	}
	if (this.neg) s = "-" + s;
	return s;
}

Bignum.prototype.normalize = function() {
	// Remove leading zeros, retaining one for exactly 0.
	//
	// This is in the inner loop of everything.
	//
	var limit = this.v.length-1;
	var v = this.v;
	for (var i = limit; i > 0; i--) if (v[i]) break;
	if (i < limit) v.length = i + 1;
	if (v.length == 1 && v[0] == 0) this.neg = false;
	if (false) {
		// debug: check that the digits are all in range
		for (var i = 0; i < v.length; i++) {
			if (v[i] < 0 || v[i] > radix) {
				alert("Bug: v[" + i + "] = " + v[i]);
				break;
			}
		}
	}
}

Bignum.prototype.absComp = function(b, e) {
	// Return (-1, 0, +1) if |this| is (<, =, >) |b * radix ^ e|
	//
	if (this.NaN || b.NaN) return Number("a"); // NaN;
	if (!e) e = 0;
	var av = this.v;
	var bv = b.v;
	var aLen = av.length;
	var bLen = bv.length;
	if (bLen == 1 && bv[0] == 0) {
		// b==0 is a special case, because shifting by e has no effect
		return (aLen > 1 || av[0] > 0 ? 1 : 0);
	}
	if (aLen < bLen + e) return -1;
	if (aLen > bLen + e) return 1;
	for (var i = aLen-1; i >= 0; i--) {
		if (i < e) {
			if (av[i] > 0) return 1;
		} else {
			if (av[i] < bv[i-e]) return -1;
			if (av[i] > bv[i-e]) return 1;
		}
	}
	return 0;
}

Bignum.prototype.comp = function(b, e) {
	// Return (-1, 0, +1) if this is (<, =, >) b * radix ^ e
	//
	if (this.neg == b.neg) {
		return (this.neg ? -this.absComp(b, e) : this.absComp(b, e));
	} else {
		if (this.NaN || b.NaN) return Number("a"); // NaN
		return (this.neg ? -1 : 1);
	}
}

Bignum.prototype.negate = function() {
	// Return negative "this"
	//
	if (this.NaN) return this;
	var res = new Bignum();
	res.v = this.v;
	res.neg = !this.neg;
	res.NaN = false;
	res.normalize();
	return res;
}

Bignum.prototype.add = function(b, e) {
	// Return "this + b * radix ^ e"
	//
	if (this.NaN) return this;
	if (b.NaN) return b;
	var av = this.v;
	var bv = b.v;
	var aLen = av.length;
	var bLen = bv.length;
	if (bLen == 1 && bv[0] == 0) return this; // negate doesn't flip 0
	if (this.neg != b.neg) return this.sub(b.negate(), e);
	if (!e) e = 0;
	var res = new Bignum();
	res.NaN = false;
	res.v = av.slice(0, aLen);
	var resV = res.v;
	var carry = 0;
	// The following loops work even if a >= bLen (which doesn't happen)
	// 1: 0..min(e,aLen) has been copied from "a"
	// 2: e..overlap is where "a" and "b" overlap, if at all
	var overlap = Math.min(aLen, e + bLen);
	for (var i = e; i < overlap; i++) {
		var n = av[i] + bv[i-e] + carry;
		if (n >= radix) {
			carry = 1;
			resV[i] = n - radix;
		} else {
			carry = 0;
			resV[i] = n;
		}
	}
	// 3: overlap..aLen is rest of "a", if any
	for (var i = overlap; i < aLen; i++) {
		var n = av[i] + carry;
		if (n >= radix) {
			carry = 1;
			resV[i] = n - radix;
		} else {
			carry = 0;
			resV[i] = n;
		}
	}
	// 3: aLen..padding is zero-fill, if any, caused by e > aLen
	var padding = Math.max(aLen, e);
	for (var i = aLen; i < padding; i++) {
		resV[i] = carry;
		carry = 0;
	}
	// 4: padding..e+bLen is rest of "b", if any
	for (var i = padding; i < e + bLen; i++) {
		var n = bv[i-e] + carry;
		if (n >= radix) {
			carry = 1;
			resV[i] = n - radix;
		} else {
			carry = 0;
			resV[i] = n;
		}
	}
	if (carry > 0) resV[res.v.length] = carry;
	res.neg = this.neg;
	res.normalize();
	return res;
}

Bignum.prototype.sub = function(b, e) {
	// Return "this - b * radix ^ e"
	//
	if (this.NaN) return this;
	if (b.NaN) return b;
	var av = this.v;
	var bv = b.v;
	var aLen = av.length;
	var bLen = bv.length;
	if (bLen == 1 && bv[0] == 0) return this; // negate doesn't flip 0
	if (this.neg != b.neg) return this.add(b.negate(), e);
	if (!e) e = 0;
	var res = new Bignum();
	res.NaN = false;
	res.v = av.slice(0, aLen);
	var resV = res.v;
	var carry = 0;
	// The following loops work even if a >= bLen (which doesn't happen)
	// For commentary, see "add".
	var overlap = Math.min(aLen, e + bLen);
	for (var i = e; i < overlap; i++) {
		var n = av[i] - bv[i-e] - carry;
		if (n < 0) {
			carry = 1;
			resV[i] = n + radix;
		} else {
			carry = 0;
			resV[i] = n;
		}
	}
	for (var i = overlap; i < aLen; i++) {
		var n = av[i] - carry;
		if (n < 0) {
			carry = 1;
			resV[i] = n + radix;
		} else {
			carry = 0;
			resV[i] = n;
		}
	}
	var padding = Math.max(aLen, e);
	for (var i = aLen; i < padding; i++) {
		resV[i] = (carry > 0 ? radix-1 : 0);
	}
	for (var i = padding; i < e + bLen; i++) {
		var n = -bv[i-e] - carry;
		if (n < 0) {
			carry = 1;
			resV[i] = n + radix;
		} else {
			carry = 0;
			resV[i] = n;
		}
	}
	res.neg = this.neg;
	if (carry > 0) {
		// Sadly, |this| < |b|
		// The answer is "-(radix^resLen) + res", so we compute
		// -(radix^resLen - res).
		// Because of "e", this is simpler than detecting the problem up
		// front and swapping the arguments.
		var resLen = res.v.length;
		carry = 0;
		for (var i = 0; i < resLen; i++) {
			var n = 0 - resV[i] - carry;
			if (n < 0) {
				carry = 1;
				resV[i] = n + radix;
			} else {
				carry = 0;
				resV[i] = n;
			}
		}
		res.neg = !res.neg;
	}
	res.normalize();
	return res;
}

function bigSlice(a, start, end) {
	// INTERNAL: return a number composed of digits [start .. end-1] of "a"
	//
	if (a.NaN) return a;
	var res = new Bignum();
	res.v = a.v.slice(start, end);
	res.NaN = false;
	res.neg = this.neg;
//	res.normalize();
//  Not needed, and surprisingly slow, with lots of calls from bigMulK
	return res;
}

function bigMulClassic(a, b) {
	// INTERNAL: return "a * b", assuming a and b are both unsigned non-NaN.
	//
	// This version is classic paper-and-pencil long multiplication,
	// with the order of the multiplications shuffled a little.
	//
	// The special case of squaring (a==b) is optimized.
	//
	var big = (a.v.length >= b.v.length ? a : b);
	var small = (big == a ? b : a);
	var bigV = big.v;
	var smallV = small.v;
	var bigLen = bigV.length;
	var smallLen = smallV.length;
	var res = new Bignum(0);
	var resV = res.v;
	for (var i = 0; i < smallLen; i++) {
		var nBig = bigV[i];
		var nSmall = smallV[i];
		var carry = 0;
		for (var j = 0; j < i; j++) {
			var n = nBig * smallV[j];
			n += resV[i+j] + (big == small ? n : nSmall * bigV[j]) + carry;
			// Note: "n" can reach 2*(radix*radix-1)
			carry = Math.floor(n / radix);
			resV[i+j] = n - carry * radix;
		}
		// Do the diagonal
		var n = nBig * nSmall + carry;
		carry = Math.floor(n / radix);
		resV[i+i] = n - carry * radix;
		resV[i+i+1] = carry;
	}
	for (var i = smallLen; i < bigLen; i++) {
		var nBig = bigV[i];
		var carry = 0;
		for (var j = 0; j < smallLen; j++) {
			var n = resV[i+j] + nBig * smallV[j] + carry;
			carry = Math.floor(n / radix);
			resV[i+j] = n - carry * radix;
		}
		resV[i+smallLen] = carry;
	}
	res.normalize();
	return res;
}

function bigMulK(a, b) {
	// INTERNAL return "a * b", assuming a and b are both unsigned non-NaN.
	//
	// This version uses the Karatsuba algorithm.
	//
	var aLen = a.v.length;
	var bLen = b.v.length;
	var half = Math.max(Math.floor(aLen/2), Math.floor(bLen/2));
	if (half < 20) {
		// For small number of digits, classic is faster
		return bigMulClassic(a, b);
	} else if (half >= aLen) {
		// "a" is too short: there's nothing to gain
		var msB = bigSlice(b, half, bLen);
		var lsB = bigSlice(b, 0, half);
		var res = bigMulK(a, lsB).add(bigMulK(a, msB), half);
		return res;
	} else if (half >= bLen) {
		// "b" is too short: there's nothing to gain
		var msA = bigSlice(a, half, aLen);
		var lsA = bigSlice(a, 0, half);
		var res = bigMulK(b, lsA).add(bigMulK(b, msA), half);
		return res;
	} else if (a == b) {
		// Karatsuba squaring
		var msA = bigSlice(a, half, aLen);
		var lsA = bigSlice(a, 0, half);
		var ms = bigMulK(msA, msA);
		var ls = bigMulK(lsA, lsA);
		var midA = msA.add(lsA);
		var mid = bigMulK(midA, midA).sub(ms).sub(ls);
		var res = ls.add(mid, half).add(ms, half+half);
		return res;
	} else {
		// Karatsuba proper
		var msA = bigSlice(a, half, aLen);
		var lsA = bigSlice(a, 0, half);
		var msB = bigSlice(b, half, bLen);
		var lsB = bigSlice(b, 0, half);
		var ms = bigMulK(msA, msB);
		var ls = bigMulK(lsA, lsB);
		var midA = msA.add(lsA);
		var midB = msB.add(lsB);
		var mid = bigMulK(midA, midB).sub(ms).sub(ls);
			// mid = (msA+lsA)*(msB+lsB)-msA*msB-lsA*lsB
			//     = msA*lsB+lsA*msB
		var res = ls.add(mid, half).add(ms, half+half);
		return res;
	}
}

Bignum.prototype.mul = function(b, m) {
	// Return "this * b" or "(this * b) mod m" depending on m.
	//
	if (this.NaN) return this;
	if (b.NaN) return b;
	if (m && m.NaN) return m;
	var res = bigMulK(this, b);
	res.neg = (this.neg != b.neg);
	res.normalize();
	return (m ? res.mod(m) : res);
}

Bignum.prototype.mulInt = function(n) {
	// Return "this * n".
	//
	// Entirely equivalent to this.mul(new Bignum(n)), but faster
	//
	var av = this.v;
	var aLen = av.length;
	if (this.NaN) return this;
	if (n > radix) {
		return this.mul(new Bignum(n));
	} else if (n < 0) {
		return this.mulInt(-n).negate();
	} else if (n == 0) {
		return bigZero;
	} else if (n == 1) {
		return this;
	} else {
		var res = new Bignum(0);
		var resV = res.v;
		var carry = 0;
		for (var i = 0; i < aLen; i++) {
			var x = n * av[i] + carry;
			carry = Math.floor(x / radix);
			resV[i] = x - carry * radix;
		}
		if (carry > 0) resV[aLen] = carry;
		// no need to normalize
		res.neg = this.neg;
		return res;
	}
}

function bigMulIntSub(a, b, n, e) {
	// INTERNAL: return "a - b * n * radix^e", unless it would be
	// negative, in which case return false.
	//
	// Assumes unsigned and non-NaN, and requires 0 <= n < radix and b != 0
	//
	// This is the inner loop, and almost all of the time, of "divmod".
	//
	var av = a.v;
	var bv = b.v;
	var aLen = av.length;
	var bLen = bv.length;
	if (n < 0 || n >= radix) {
		return new Bignum("internal error in bigMulIntSub");
	} else if (n == 0) {
		return a;
	} else if (e + bLen > aLen) {
		return false; // clearly too big.  This check is relied on later.
	} else {
		var res = new Bignum();
		res.v = av.slice(0, e);
		res.NaN = false;
		var resV = res.v;
		var carry = 0;
		for (var i = e; i < e + bLen; i++) {
			var x = av[i] - n * bv[i-e] - carry;
			if (x < 0) {
				carry = Math.ceil((-x) / radix);
				resV[i] = x + carry * radix;
			} else {
				carry = 0;
				resV[i] = x;
			}
		}
		for (var i = e + bLen; i < aLen; i++) {
			var x = av[i] - carry;
			if (x < 0) {
				carry = 1;
				resV[i] = x + radix;
			} else {
				carry = 0;
				resV[i] = x;
			}
		}
		// The subtractand is too big iff we still have a carry.
		if (carry > 0) return false;
		res.normalize();
		return res;
	}
}

var divModIter = 0;  // debug

Bignum.prototype.divmod = function(b, mod) {
	// Return "this / b" or "this mod b", depending on "mod".
	// "this / b" is rounded toward zero.
	// "this mod b" equals "this - (this / b) * b"
	//
	// We use classic shift/test/subtract long division.  This is just like
	// it's taught in grade school, but our number system is base "radix"
	// instead of base 10.
	//
	// We take advantage of the "e" argument of comp and sub to avoid
	// actually performing the shifts for each digit.
	//
	if (this.NaN) return this;
	if (b.NaN) return b;
	if (b.absComp(bigZero) == 0) return new Bignum("divide by zero");
	var res = (mod ? null : new Bignum(0));
	var av = this.v;
	var bv = b.v;
	var aLen = av.length;
	var bLen = bv.length;
	var msDigits = bv[bLen-1] * radix + (bLen-2 >= 0 ? bv[bLen-2] : 0);
	// Do the work in positive, then fix up signs at the end
	var rem = (this.neg ? this.negate() : this);
	var divisor = (b.neg ? b.negate() : b);
	for (var e = aLen - bLen; e >= 0; e--) {
		var remV = rem.v;
		var remLen = remV.length;
		var thisDigit;
		var divDigit0 = (remLen > bLen+e ? remV[bLen+e] : 0);
		var divDigit1 = (remLen > bLen+e-1 && bLen+e-1 >= 0 ?
												remV[bLen+e-1] : 0);
		var divDigit2 = (remLen > bLen+e-2 && bLen+e-2 >= 0 ?
												remV[bLen+e-2] : 0);
		var divDigits = (divDigit0 * radix + divDigit1) * radix + divDigit2;
		var thisDigit = Math.floor(divDigits / msDigits);
			// This is just an estimate, though almost always correct.
			// If it's wrong it's almost always high; but when divDigits
			// (which is a float) exceeds maxint, rounding can make
			// thisDigit one too low.
		while (true) {
			var nextRem = bigMulIntSub(rem, divisor, thisDigit, e);
			if (nextRem) { rem = nextRem; break; }
			thisDigit--;
		}
		for (var incrs = 0; true; incrs++) {
			if (rem.absComp(divisor, e) < 0) break; // rem < divisor, good.
			// "thisDigit" was too low
			rem = rem.sub(divisor, e);
			thisDigit++;
			if (incrs > 1) {
				alert("Bug: under-estimated digit in divmod: " + thisDigit);
				return new Bignum("internal error");
			}
		}
		if (!mod) res.v[e] = thisDigit;
	}
	if (rem.absComp(divisor) >= 0) {
		alert("Bug: remainder is too large");
		return new Bignum("internal error");
	}
	// Fix up signs, and normalize "res".
	if (!mod) {
		res.neg = (this.neg != b.neg);
		res.normalize();
	}
	if (this.neg) rem = rem.negate();
	return (mod ? rem : res);
}

Bignum.prototype.div = function(b) {
	// Return "this / b"
	//
	return this.divmod(b, false);
}

Bignum.prototype.mod = function(b) {
	// Return "this mod b"
	//
	return this.divmod(b, true);
}

function bigPowLS(a, b, m) {
	// INTERNAL: return "a^b mod m", assuming b is unsigned non-Nan and m
	// is non-NaN.
	//
	// This variant does l.s. bit first, with no assumption about "radix"
	//
	var exp = a;
	var res = bigOne;
	var t = bigClone(b);
	while (t.length > 0) {
		var limit = t.length-1;
		if (t[limit] == 0) {
			t.length = limit;
		} else {
			var rem = bigDivInt(t, 2);
			if (rem > 0) {
				res = res.mul(exp, m);
				if (res.v.length * radixBits > 500000) return new Bignum(
						"infeasible, too large");
			}
			exp = exp.mul(exp, m);
		}
	}
	return res;
}

function bigPowMS(a, b, m) {
	// INTERNAL: return "a^b mod m", assuming b is unsigned non-Nan and m
	// is non-NaN.
	//
	// This variant does m.s. bit first, and assumes "radix" is a power of
	// 2.  It should be faster for small "a" and a wash otherwise.
	//
	var res = bigOne;
	for (var i = b.v.length-1; i >= 0; i--) {
		var j = radix;
		var digit = b.v[i];
		while ((j = Math.floor(j / 2)) > 0) {
			res = res.mul(res, m);
			if (j <= digit) {
				digit -= j;
				res = res.mul(a, m);
			}
			if (res.NaN) return res;
			if (res.v.length * radixBits > 500000) return new Bignum(
					"infeasible, too large");
		}
	}		
	return res;
}

Bignum.prototype.pow = function(b, m) {
	// Return "this ^ b" or "(this ^ b) mod m" depending on m.
	// Note that, following EcmaScript, NaN.pow(0) == 1.
	//
	if (b.NaN) return b;
	if (m && m.NaN) return m;
	if (b.neg) return new Bignum("negative exponent");
	return bigPowMS(this, b, m);
}

Bignum.prototype.gcd = function(b) {
	// Return GCD(this, b).
	//
	// A trivial application of Euclid's algorithm.
	//
	if (this.NaN) return this;
	if (b.NaN) return b;
	if (this.neg) return this.negate().gcd(b);
	if (b.neg) return this.gcd(b.negate());
	var big = (this.absComp(b) > 0 ? this : b);
	var small = (big == this ? b : this);
	while (small.absComp(bigZero) != 0) {
		var rem = big.mod(small);
		big = small;
		small = rem;
	}
	return big;
}

Bignum.prototype.inverse = function(b) {
	// Return "c" s.t. "(this * c) mod b == 1", i.e., the modular inverse.
	//
	// A special case of the extended Euclidean algorithm.
	// Requires GCD(this, b) == 1.
	//
	if (this.NaN) return this;
	if (b.NaN) return b;
	if (this.neg || b.neg)
		return new Bignum("modular inverse of negative number");
	var x = bigOne;
	var lastX = bigZero;
	var big = b;
	var small = (this.absComp(b) < 0 ? this : this.mod(b));
	while (small.absComp(bigZero) != 0) {
		var quotient = big.div(small);
		var rem = big.mod(small);
		big = small;
		small = rem;
		var t = x;
		x = lastX.sub(quotient.mul(x));
		lastX = t;
	}
	if (big.absComp(bigOne) != 0)
		return new Bignum("modular inverse with GCD not 1");
	if (lastX.neg) lastX = lastX.add(b);
	return lastX;
}

Bignum.prototype.rabin = function(nMinus1, s, d, a) {
	// Inner loop of Miller-Rabin primality test.
	// Return true if "this" appears to be prime.
	//
	// First, for convenience, ignore a <= 1 and a >= n-1
	if (a.absComp(bigOne) <= 0 || a.absComp(nMinus1) >= 0) return true;
	var p = a.pow(d, this);
	if (p.absComp(bigOne) == 0) return true;
	if (p.absComp(nMinus1) == 0) return true;
	for (var i = 0; i < s-1; i++) {
		p = p.mul(p, this);
		if (p.absComp(bigOne) == 0) return false;
		if (p.absComp(nMinus1) == 0) return true;
	}
	return false;
}

var jaeschke = new Bignum("3825123056546413051");
	// Miller-Rabin is deterministic using first 9 primes below this
	// http://oeis.org/A014233
	// G. Jaeschke, Mathematics of Computation, 61 (1993), 915-926
var jaeschkeN = 9;

var smallPrimes = null; // A cache, filled in "isPrime"
var jPrimes = null;     // Similar, with Bignum values, for Jaeschke test

function eratosthenes(n) {
	// Return an array containing the primes less than n
	//
	// The sieve of Eratosthenes
	//
	var res = new Array();
	var sieve = new Array();
	var p = 2;
	while (p * p < n) {
		res[res.length] = p;
		for (var multiple = p; multiple < n; multiple += p) {
			sieve[multiple] = true;
		}
		while (sieve[p]) p++;
	}
	for (; p < n; p++) if (!sieve[p]) res[res.length] = p;
	return res;
}

Bignum.prototype.isPrime = function() {
	// Return true iff "this" is probably prime, as indicated by running
	// the Miller-Rabin probabalistic primality test for confidence 1:2^80
	//
	// For numbers < "jaeschke", the result is true, not probabalistic.
	//
	if (this.NaN) return false;
	if (this.absComp(bigOne) <= 0) return false;

	if (!smallPrimes) {
		smallPrimes = eratosthenes(100 * 1000);
		jPrimes = new Array();
		for (var p = 0; p < jaeschkeN; p++) {
			jPrimes[p] = new Bignum(smallPrimes[p]);
		}
	}
	for (var p = 0; p < smallPrimes.length; p++) {
		var smallP = smallPrimes[p];
		if (this.v.length == 1 && this.v[0] < smallP * smallP) return true;
		// Do a check divide using schoolroom short division;
		// assumes smallP < radix, roughly
		var carry = 0;
		for (var i = this.v.length - 1; i >= 0; i--) {
			var dividend = (carry * radix + this.v[i]);
			var quotient = Math.floor(dividend / smallP);
			carry = dividend - quotient * smallP;
		}
		if (carry == 0) return false;
	}

	// First, setup n-1, s and d, s.t. d is odd and d*2^s = n-1
	var nMinus1 = this.sub(bigOne);
	var s = 0;
	var t = bigClone(nMinus1);
	while (true) { // note that nMinus1 > 0, so we always terminate.
		var limit = t.length-1;
		if (t[limit] == 0) {
			t.length = limit;
		} else if (bigDivInt(t, 2) > 0) {
			break;
		} else {
			s++;
		}
	}
	var d = new Bignum();
	d.v = t;
	d.NaN = false;
	d.normalize();
	// We overshot: we really want d = d*2+1
	d = d.mulInt(2).add(bigOne);
	// Done creating nMinus1, s, and d

	if (this.absComp(jaeschke) < 0) {
		for (var p = 0; p < jaeschkeN; p++) {
			if (!this.rabin(nMinus1, s, d, jPrimes[p])) return false;
		}
	} else {
		var thisBits = (this.v.length - 1) * radixBits; // lower bound
		// We choose "trials" to achieve an error probability 1 in 2^80,
		// using tables 4.3 and 4.4 of Handbook of Applied Cryptography.
		var trials = (thisBits >= 600 ? 2 :  // first, ones from table 4.3
					  thisBits >= 550 ? 4 :
					  thisBits >= 500 ? 5 :
					  thisBits >= 400 ? 6 :
					  thisBits >= 350 ? 7 :
					  thisBits >= 300 ? 9 :
					  thisBits >= 250 ? 12 : // then, ones from table 4.4
					  thisBits >= 200 ? 15 :
					  thisBits >= 150 ? 18 :
					  thisBits >= 100 ? 27 :
					  40);                   // finally, the basic (1/4)^t
		for (var i = 0; i < trials; i++) {
			var a = nMinus1.randomLT();
			if (!this.rabin(nMinus1, s, d, a)) return false;
			if (a.absComp(bigOne) <= 0) i--; // it didn't count
		}
	}
	return true;
}

Bignum.prototype.randomLT = function() {
	// Returns a fairly random positive Bignum less than |this|
	//
	if (this.absComp(bigZero) == 0) return new Bignum("random with 0 bits");
	var res = new Bignum(0);
	var max = this.sub(bigOne);
	var lt = false; // true iff some more significant digit is < max
	for (var i = max.v.length-1; i >= 0; i--) {
		var digit;
		if (lt) {
			while ((digit = Math.floor(Math.random() * radix)) >= radix) { }
		} else {
			while((digit = Math.floor(Math.random() * (max.v[i] + 1))) >
							max.v[i]) { } // loop for rounding errors
			if (digit < max.v[i]) lt = true;
		}
		res.v[i] = digit;
	}
	res.normalize();
	return res;
}

function bigRandom(bits) {
	// Returns a fairly random Bignum uniformly distributed in [0..2^bits).
	//
	if (!bits || bits <= 0) return new Bignum("random with 0 bits");
	return bigTwo.pow(new Bignum(bits)).randomLT();
}

function bigPrime(bits) {
	// Returns a fairly random probable prime uniformly distributed in
	// [2^(bits-1)..2^bits)
	//
	// Note that, unlike "bigRandom", the top bit is always 1.
	//
	if (!bits || bits <= 0) return new Bignum("prime with 0 bits");
	var top = bigTwo.pow(new Bignum(bits-1));
	for (var i = 0; i < 5000; i++) {
		var r = top.randomLT().add(top);
		if (r.isPrime()) return r;
	}
	return new Bignum("no prime found");
}