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

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


// 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.
//
// In general, our semantics mirror the fine print of EcmaScript Number
// semantics (except for the omission of infinity and negative zero):
//   errors are convert into "NaN";
//   "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.
//
// Overall, the performance is moderately good, and at least as good as any
// other Javascript implementation that I've seen.  Multiply uses the
// Karatsuba algorithm for numbers large enough for it to be faster.  Divmod
// uses classic schoolroom long division.
//
// The FFT multiplication algorithm is asymptotically faster than Karatsuba,
// but it wouldn't help until well beyond 1000-bit numbers.  Neither Newton-
// Raphson division nor Montgomery reduction (for a**b mod m) help unless
// you have much faster multiply.  The primary candidate for fine-grain
// optimization is the inner loop of bigMulIntSub.
//
// The public entry points are the constructor "Bignum", the method calls
// of a Bignum object, and the functions "bigRandom" and "bigPrime".  The
// other functions are internal to the implementation.
//
// This library is provided, if at all, entirely "as is"; use it at your
// own risk.  Don't even dream of using the prime generator for serious
// cryptographic keys (the primes aren't random enough, apart from anything
// else).


var radix = 1024*1024;
	// EcmaScript represents numbers as 64-bit IEEE floats, with exact
	// integers up to 2**53 (roughly 9 * 10**15).
	//
	// That should allow radix 2**26; however this fails under Safari
	// 4.0.4 on MacOS 10.5 on PowerPC.  2**20 seems to work correctly.
	//
	// As coded below, we require exact representation of integers up to
	// 2*(radix*radix-1) ... see the inner loop of "bigMulClassic".
	//
	// "bigPowMS" assumes radix is a power of 2.
	//
var radixBits = 24;
	// 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 = true;
		} else {
			this.v = temp.v;
		}
	}
}

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

Bignum.prototype.clone = function() {
	// Returns a copy of this number's integer array, primarily for use by
	// bigDivInt.
	//
	return this.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 Number("a").toString();
	var clump = 10000000;
	var digits = 7;
	var t = this.clone();
	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;
}

Bignum.prototype.slice = function(start, end) {
	// Return a number composed of digits [start .. end-1] of "this"
	//
	if (this.NaN) return this;
	var res = new Bignum();
	res.v = this.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.
	//
	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 < 12) {
		// 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 = b.slice(half, bLen);
		var lsB = b.slice(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 = a.slice(half, aLen);
		var lsA = a.slice(0, half);
		var res = bigMulK(b, lsA).add(bigMulK(b, msA), half);
		return res;
	} else if (a == b) {
		// Karatsuba squaring
		var msA = a.slice(half, aLen);
		var lsA = a.slice(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 = a.slice(half, aLen);
		var lsA = a.slice(0, half);
		var msB = b.slice(half, bLen);
		var lsB = b.slice(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 {
		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(); // caller error
	} 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();
	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();
			}
		}
		if (!mod) res.v[e] = thisDigit;
	}
	if (rem.absComp(divisor) >= 0) {
		alert("Bug: remainder is too large");
		return new Bignum();
	}
	// 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 = b.clone();
	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);
			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);
			}
		}
	}		
	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();
	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();
	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();
	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(341550071728321); // 3.4 * 10*14
	// Miller-Rabin is deterministic using first 7 primes below this
var jaeschkeN = 7;
var smallPrimesRaw = new Array( // primes below 1000
	  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);
var smallPrimes = new Array(); // A cache, filled in "isPrime"

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
	//
	if (this.NaN) return false;

	// First, the easy ones ...
	if (this.absComp(bigOne) <= 0) return false;
	for (var p = 0; p < smallPrimesRaw.length; p++) {
		var pNum = smallPrimes[p];
		if (!pNum) pNum = (smallPrimes[p] = new Bignum(smallPrimesRaw[p]));
		if (this.absComp(pNum) <= 0) return true;
		if (this.mod(pNum).absComp(bigZero) == 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 = nMinus1.clone();
	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

	// It's always worth trying a==2
	if (!this.rabin(nMinus1, s, d, smallPrimes[0])) return false;

	if (this.absComp(jaeschke) < 0) {
		for (var p = 1; p < jaeschkeN; p++) {
			if (!this.rabin(nMinus1, s, d, smallPrimes[p])) return false;
		}
	} else {
		var thisBits = this.v.length * radixBits; // high guess
		// We choose "trials" to achieve an error probability 1 in 2**80,
		// using numbers from table 4.4 of Handbook of Applied Cryptography.
		trials = (thisBits < 150 ? 27 : thisBits < 250 ? 15 :
			thisBits < 550 ? 6 : thisBits < 1300 ? 3 : 2);
		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|
	//
	var res = new Bignum(0);
	for (var i = 0; i < this.v.length-1; i++) {
		res.v[i] = Math.floor(Math.random() * radix);
	}
	res.v[this.v.length-1] =
					Math.floor(Math.random() * this.v[this.v.length-1]);
	res.normalize();
	return res;
}

function bigRandom(thisBits) {
	// Returns a fairly random Bignum having approximately thisBits bits.
	//
	if (!thisBits || thisBits <= 0) return new Bignum();
	var res = new Bignum(0);
	var chunkLen = 25;
	var chunkMax = 1024 * 1024 * 32;
	var remaining = thisBits;
	while (remaining >= chunkLen) {
		res = res.mulInt(chunkMax);
		res = res.add(new Bignum(Math.floor(Math.random() * chunkMax)));
		remaining -= chunkLen;
	}
	var tail = 1;
	while (remaining > 0) { tail = tail * 2; remaining--; }
	res = res.mulInt(tail);
	res = res.add(new Bignum(Math.floor(Math.random() * tail)));
	return res;
}

function bigPrime(thisBits) {
	// Returns a fairly random probable prime with about n bits
	//
	if (!thisBits || thisBits <= 0) return new Bignum();
	for (var i = 0; i < 1000; i++) {
		var r = bigRandom(thisBits);
		if (r.isPrime()) return r;
	}
	return new Bignum();
}