begin crypt

18 files changed, 3162 insertions, 0 deletions

diff --git a/bigint/BigInteger.cc b/bigint/BigInteger.cc
new file mode 100644
index 0000000..3b23aa1
--- /dev/null
+++ b/bigint/BigInteger.cc
@@ -0,0 +1,405 @@
+#include "BigInteger.hh"
+
+void BigInteger::operator =(const BigInteger &x) {
+	// Calls like a = a have no effect
+	if (this == &x)
+		return;
+	// Copy sign
+	sign = x.sign;
+	// Copy the rest
+	mag = x.mag;
+}
+
+BigInteger::BigInteger(const Blk *b, Index blen, Sign s) : mag(b, blen) {
+	switch (s) {
+	case zero:
+		if (!mag.isZero())
+			throw "BigInteger::BigInteger(const Blk *, Index, Sign): Cannot use a sign of zero with a nonzero magnitude";
+		sign = zero;
+		break;
+	case positive:
+	case negative:
+		// If the magnitude is zero, force the sign to zero.
+		sign = mag.isZero() ? zero : s;
+		break;
+	default:
+		/* g++ seems to be optimizing out this case on the assumption
+		 * that the sign is a valid member of the enumeration.  Oh well. */
+		throw "BigInteger::BigInteger(const Blk *, Index, Sign): Invalid sign";
+	}
+}
+
+BigInteger::BigInteger(const BigUnsigned &x, Sign s) : mag(x) {
+	switch (s) {
+	case zero:
+		if (!mag.isZero())
+			throw "BigInteger::BigInteger(const BigUnsigned &, Sign): Cannot use a sign of zero with a nonzero magnitude";
+		sign = zero;
+		break;
+	case positive:
+	case negative:
+		// If the magnitude is zero, force the sign to zero.
+		sign = mag.isZero() ? zero : s;
+		break;
+	default:
+		/* g++ seems to be optimizing out this case on the assumption
+		 * that the sign is a valid member of the enumeration.  Oh well. */
+		throw "BigInteger::BigInteger(const BigUnsigned &, Sign): Invalid sign";
+	}
+}
+
+/* CONSTRUCTION FROM PRIMITIVE INTEGERS
+ * Same idea as in BigUnsigned.cc, except that negative input results in a
+ * negative BigInteger instead of an exception. */
+
+// Done longhand to let us use initialization.
+BigInteger::BigInteger(unsigned long  x) : mag(x) { sign = mag.isZero() ? zero : positive; }
+BigInteger::BigInteger(unsigned int   x) : mag(x) { sign = mag.isZero() ? zero : positive; }
+BigInteger::BigInteger(unsigned short x) : mag(x) { sign = mag.isZero() ? zero : positive; }
+
+// For signed input, determine the desired magnitude and sign separately.
+
+namespace {
+	template <class X, class UX>
+	BigInteger::Blk magOf(X x) {
+		/* UX(...) cast needed to stop short(-2^15), which negates to
+		 * itself, from sign-extending in the conversion to Blk. */
+		return BigInteger::Blk(x < 0 ? UX(-x) : x);
+	}
+	template <class X>
+	BigInteger::Sign signOf(X x) {
+		return (x == 0) ? BigInteger::zero
+			: (x > 0) ? BigInteger::positive
+			: BigInteger::negative;
+	}
+}
+
+BigInteger::BigInteger(long  x) : sign(signOf(x)), mag(magOf<long , unsigned long >(x)) {}
+BigInteger::BigInteger(int   x) : sign(signOf(x)), mag(magOf<int  , unsigned int  >(x)) {}
+BigInteger::BigInteger(short x) : sign(signOf(x)), mag(magOf<short, unsigned short>(x)) {}
+
+// CONVERSION TO PRIMITIVE INTEGERS
+
+/* Reuse BigUnsigned's conversion to an unsigned primitive integer.
+ * The friend is a separate function rather than
+ * BigInteger::convertToUnsignedPrimitive to avoid requiring BigUnsigned to
+ * declare BigInteger. */
+template <class X>
+inline X convertBigUnsignedToPrimitiveAccess(const BigUnsigned &a) {
+	return a.convertToPrimitive<X>();
+}
+
+template <class X>
+X BigInteger::convertToUnsignedPrimitive() const {
+	if (sign == negative)
+		throw "BigInteger::to<Primitive>: "
+			"Cannot convert a negative integer to an unsigned type";
+	else
+		return convertBigUnsignedToPrimitiveAccess<X>(mag);
+}
+
+/* Similar to BigUnsigned::convertToPrimitive, but split into two cases for
+ * nonnegative and negative numbers. */
+template <class X, class UX>
+X BigInteger::convertToSignedPrimitive() const {
+	if (sign == zero)
+		return 0;
+	else if (mag.getLength() == 1) {
+		// The single block might fit in an X.  Try the conversion.
+		Blk b = mag.getBlock(0);
+		if (sign == positive) {
+			X x = X(b);
+			if (x >= 0 && Blk(x) == b)
+				return x;
+		} else {
+			X x = -X(b);
+			/* UX(...) needed to avoid rejecting conversion of
+			 * -2^15 to a short. */
+			if (x < 0 && Blk(UX(-x)) == b)
+				return x;
+		}
+		// Otherwise fall through.
+	}
+	throw "BigInteger::to<Primitive>: "
+		"Value is too big to fit in the requested type";
+}
+
+unsigned long  BigInteger::toUnsignedLong () const { return convertToUnsignedPrimitive<unsigned long >       (); }
+unsigned int   BigInteger::toUnsignedInt  () const { return convertToUnsignedPrimitive<unsigned int  >       (); }
+unsigned short BigInteger::toUnsignedShort() const { return convertToUnsignedPrimitive<unsigned short>       (); }
+long           BigInteger::toLong         () const { return convertToSignedPrimitive  <long , unsigned long> (); }
+int            BigInteger::toInt          () const { return convertToSignedPrimitive  <int  , unsigned int>  (); }
+short          BigInteger::toShort        () const { return convertToSignedPrimitive  <short, unsigned short>(); }
+
+// COMPARISON
+BigInteger::CmpRes BigInteger::compareTo(const BigInteger &x) const {
+	// A greater sign implies a greater number
+	if (sign < x.sign)
+		return less;
+	else if (sign > x.sign)
+		return greater;
+	else switch (sign) {
+		// If the signs are the same...
+	case zero:
+		return equal; // Two zeros are equal
+	case positive:
+		// Compare the magnitudes
+		return mag.compareTo(x.mag);
+	case negative:
+		// Compare the magnitudes, but return the opposite result
+		return CmpRes(-mag.compareTo(x.mag));
+	default:
+		throw "BigInteger internal error";
+	}
+}
+
+/* COPY-LESS OPERATIONS
+ * These do some messing around to determine the sign of the result,
+ * then call one of BigUnsigned's copy-less operations. */
+
+// See remarks about aliased calls in BigUnsigned.cc .
+#define DTRT_ALIASED(cond, op) \
+	if (cond) { \
+		BigInteger tmpThis; \
+		tmpThis.op; \
+		*this = tmpThis; \
+		return; \
+	}
+
+void BigInteger::add(const BigInteger &a, const BigInteger &b) {
+	DTRT_ALIASED(this == &a || this == &b, add(a, b));
+	// If one argument is zero, copy the other.
+	if (a.sign == zero)
+		operator =(b);
+	else if (b.sign == zero)
+		operator =(a);
+	// If the arguments have the same sign, take the
+	// common sign and add their magnitudes.
+	else if (a.sign == b.sign) {
+		sign = a.sign;
+		mag.add(a.mag, b.mag);
+	} else {
+		// Otherwise, their magnitudes must be compared.
+		switch (a.mag.compareTo(b.mag)) {
+		case equal:
+			// If their magnitudes are the same, copy zero.
+			mag = 0;
+			sign = zero;
+			break;
+			// Otherwise, take the sign of the greater, and subtract
+			// the lesser magnitude from the greater magnitude.
+		case greater:
+			sign = a.sign;
+			mag.subtract(a.mag, b.mag);
+			break;
+		case less:
+			sign = b.sign;
+			mag.subtract(b.mag, a.mag);
+			break;
+		}
+	}
+}
+
+void BigInteger::subtract(const BigInteger &a, const BigInteger &b) {
+	// Notice that this routine is identical to BigInteger::add,
+	// if one replaces b.sign by its opposite.
+	DTRT_ALIASED(this == &a || this == &b, subtract(a, b));
+	// If a is zero, copy b and flip its sign.  If b is zero, copy a.
+	if (a.sign == zero) {
+		mag = b.mag;
+		// Take the negative of _b_'s, sign, not ours.
+		// Bug pointed out by Sam Larkin on 2005.03.30.
+		sign = Sign(-b.sign);
+	} else if (b.sign == zero)
+		operator =(a);
+	// If their signs differ, take a.sign and add the magnitudes.
+	else if (a.sign != b.sign) {
+		sign = a.sign;
+		mag.add(a.mag, b.mag);
+	} else {
+		// Otherwise, their magnitudes must be compared.
+		switch (a.mag.compareTo(b.mag)) {
+			// If their magnitudes are the same, copy zero.
+		case equal:
+			mag = 0;
+			sign = zero;
+			break;
+			// If a's magnitude is greater, take a.sign and
+			// subtract a from b.
+		case greater:
+			sign = a.sign;
+			mag.subtract(a.mag, b.mag);
+			break;
+			// If b's magnitude is greater, take the opposite
+			// of b.sign and subtract b from a.
+		case less:
+			sign = Sign(-b.sign);
+			mag.subtract(b.mag, a.mag);
+			break;
+		}
+	}
+}
+
+void BigInteger::multiply(const BigInteger &a, const BigInteger &b) {
+	DTRT_ALIASED(this == &a || this == &b, multiply(a, b));
+	// If one object is zero, copy zero and return.
+	if (a.sign == zero || b.sign == zero) {
+		sign = zero;
+		mag = 0;
+		return;
+	}
+	// If the signs of the arguments are the same, the result
+	// is positive, otherwise it is negative.
+	sign = (a.sign == b.sign) ? positive : negative;
+	// Multiply the magnitudes.
+	mag.multiply(a.mag, b.mag);
+}
+
+/*
+ * DIVISION WITH REMAINDER
+ * Please read the comments before the definition of
+ * `BigUnsigned::divideWithRemainder' in `BigUnsigned.cc' for lots of
+ * information you should know before reading this function.
+ *
+ * Following Knuth, I decree that x / y is to be
+ * 0 if y==0 and floor(real-number x / y) if y!=0.
+ * Then x % y shall be x - y*(integer x / y).
+ *
+ * Note that x = y * (x / y) + (x % y) always holds.
+ * In addition, (x % y) is from 0 to y - 1 if y > 0,
+ * and from -(|y| - 1) to 0 if y < 0.  (x % y) = x if y = 0.
+ *
+ * Examples: (q = a / b, r = a % b)
+ *	a	b	q	r
+ *	===	===	===	===
+ *	4	3	1	1
+ *	-4	3	-2	2
+ *	4	-3	-2	-2
+ *	-4	-3	1	-1
+ */
+void BigInteger::divideWithRemainder(const BigInteger &b, BigInteger &q) {
+	// Defend against aliased calls;
+	// same idea as in BigUnsigned::divideWithRemainder .
+	if (this == &q)
+		throw "BigInteger::divideWithRemainder: Cannot write quotient and remainder into the same variable";
+	if (this == &b || &q == &b) {
+		BigInteger tmpB(b);
+		divideWithRemainder(tmpB, q);
+		return;
+	}
+
+	// Division by zero gives quotient 0 and remainder *this
+	if (b.sign == zero) {
+		q.mag = 0;
+		q.sign = zero;
+		return;
+	}
+	// 0 / b gives quotient 0 and remainder 0
+	if (sign == zero) {
+		q.mag = 0;
+		q.sign = zero;
+		return;
+	}
+
+	// Here *this != 0, b != 0.
+
+	// Do the operands have the same sign?
+	if (sign == b.sign) {
+		// Yes: easy case.  Quotient is zero or positive.
+		q.sign = positive;
+	} else {
+		// No: harder case.  Quotient is negative.
+		q.sign = negative;
+		// Decrease the magnitude of the dividend by one.
+		mag--;
+		/*
+		 * We tinker with the dividend before and with the
+		 * quotient and remainder after so that the result
+		 * comes out right.  To see why it works, consider the following
+		 * list of examples, where A is the magnitude-decreased
+		 * a, Q and R are the results of BigUnsigned division
+		 * with remainder on A and |b|, and q and r are the
+		 * final results we want:
+		 *
+		 *	a	A	b	Q	R	q	r
+		 *	-3	-2	3	0	2	-1	0
+		 *	-4	-3	3	1	0	-2	2
+		 *	-5	-4	3	1	1	-2	1
+		 *	-6	-5	3	1	2	-2	0
+		 *
+		 * It appears that we need a total of 3 corrections:
+		 * Decrease the magnitude of a to get A.  Increase the
+		 * magnitude of Q to get q (and make it negative).
+		 * Find r = (b - 1) - R and give it the desired sign.
+		 */
+	}
+
+	// Divide the magnitudes.
+	mag.divideWithRemainder(b.mag, q.mag);
+
+	if (sign != b.sign) {
+		// More for the harder case (as described):
+		// Increase the magnitude of the quotient by one.
+		q.mag++;
+		// Modify the remainder.
+		mag.subtract(b.mag, mag);
+		mag--;
+	}
+
+	// Sign of the remainder is always the sign of the divisor b.
+	sign = b.sign;
+
+	// Set signs to zero as necessary.  (Thanks David Allen!)
+	if (mag.isZero())
+		sign = zero;
+	if (q.mag.isZero())
+		q.sign = zero;
+
+	// WHEW!!!
+}
+
+// Negation
+void BigInteger::negate(const BigInteger &a) {
+	DTRT_ALIASED(this == &a, negate(a));
+	// Copy a's magnitude
+	mag = a.mag;
+	// Copy the opposite of a.sign
+	sign = Sign(-a.sign);
+}
+
+// INCREMENT/DECREMENT OPERATORS
+
+// Prefix increment
+void BigInteger::operator ++() {
+	if (sign == negative) {
+		mag--;
+		if (mag == 0)
+			sign = zero;
+	} else {
+		mag++;
+		sign = positive; // if not already
+	}
+}
+
+// Postfix increment: same as prefix
+void BigInteger::operator ++(int) {
+	operator ++();
+}
+
+// Prefix decrement
+void BigInteger::operator --() {
+	if (sign == positive) {
+		mag--;
+		if (mag == 0)
+			sign = zero;
+	} else {
+		mag++;
+		sign = negative;
+	}
+}
+
+// Postfix decrement: same as prefix
+void BigInteger::operator --(int) {
+	operator --();
+}
+
diff --git a/bigint/BigInteger.hh b/bigint/BigInteger.hh
new file mode 100644
index 0000000..cf6e910
--- /dev/null
+++ b/bigint/BigInteger.hh
@@ -0,0 +1,215 @@
+#ifndef BIGINTEGER_H
+#define BIGINTEGER_H
+
+#include "BigUnsigned.hh"
+
+/* A BigInteger object represents a signed integer of size limited only by
+ * available memory.  BigUnsigneds support most mathematical operators and can
+ * be converted to and from most primitive integer types.
+ *
+ * A BigInteger is just an aggregate of a BigUnsigned and a sign.  (It is no
+ * longer derived from BigUnsigned because that led to harmful implicit
+ * conversions.) */
+class BigInteger {
+
+public:
+	typedef BigUnsigned::Blk Blk;
+	typedef BigUnsigned::Index Index;
+	typedef BigUnsigned::CmpRes CmpRes;
+	static const CmpRes
+		less    = BigUnsigned::less   ,
+		equal   = BigUnsigned::equal  ,
+		greater = BigUnsigned::greater;
+	// Enumeration for the sign of a BigInteger.
+	enum Sign { negative = -1, zero = 0, positive = 1 };
+
+protected:
+	Sign sign;
+	BigUnsigned mag;
+
+public:
+	// Constructs zero.
+	BigInteger() : sign(zero), mag() {}
+
+	// Copy constructor
+	BigInteger(const BigInteger &x) : sign(x.sign), mag(x.mag) {};
+
+	// Assignment operator
+	void operator=(const BigInteger &x);
+
+	// Constructor that copies from a given array of blocks with a sign.
+	BigInteger(const Blk *b, Index blen, Sign s);
+
+	// Nonnegative constructor that copies from a given array of blocks.
+	BigInteger(const Blk *b, Index blen) : mag(b, blen) {
+		sign = mag.isZero() ? zero : positive;
+	}
+
+	// Constructor from a BigUnsigned and a sign
+	BigInteger(const BigUnsigned &x, Sign s);
+
+	// Nonnegative constructor from a BigUnsigned
+	BigInteger(const BigUnsigned &x) : mag(x) {
+		sign = mag.isZero() ? zero : positive;
+	}
+
+	// Constructors from primitive integer types
+	BigInteger(unsigned long  x);
+	BigInteger(         long  x);
+	BigInteger(unsigned int   x);
+	BigInteger(         int   x);
+	BigInteger(unsigned short x);
+	BigInteger(         short x);
+
+	/* Converters to primitive integer types
+	 * The implicit conversion operators caused trouble, so these are now
+	 * named. */
+	unsigned long  toUnsignedLong () const;
+	long           toLong         () const;
+	unsigned int   toUnsignedInt  () const;
+	int            toInt          () const;
+	unsigned short toUnsignedShort() const;
+	short          toShort        () const;
+protected:
+	// Helper
+	template <class X> X convertToUnsignedPrimitive() const;
+	template <class X, class UX> X convertToSignedPrimitive() const;
+public:
+
+	// ACCESSORS
+	Sign getSign() const { return sign; }
+	/* The client can't do any harm by holding a read-only reference to the
+	 * magnitude. */
+	const BigUnsigned &getMagnitude() const { return mag; }
+
+	// Some accessors that go through to the magnitude
+	Index getLength() const { return mag.getLength(); }
+	Index getCapacity() const { return mag.getCapacity(); }
+	Blk getBlock(Index i) const { return mag.getBlock(i); }
+	bool isZero() const { return sign == zero; } // A bit special
+
+	// COMPARISONS
+
+	// Compares this to x like Perl's <=>
+	CmpRes compareTo(const BigInteger &x) const;
+
+	// Ordinary comparison operators
+	bool operator ==(const BigInteger &x) const {
+		return sign == x.sign && mag == x.mag;
+	}
+	bool operator !=(const BigInteger &x) const { return !operator ==(x); };
+	bool operator < (const BigInteger &x) const { return compareTo(x) == less   ; }
+	bool operator <=(const BigInteger &x) const { return compareTo(x) != greater; }
+	bool operator >=(const BigInteger &x) const { return compareTo(x) != less   ; }
+	bool operator > (const BigInteger &x) const { return compareTo(x) == greater; }
+
+	// OPERATORS -- See the discussion in BigUnsigned.hh.
+	void add     (const BigInteger &a, const BigInteger &b);
+	void subtract(const BigInteger &a, const BigInteger &b);
+	void multiply(const BigInteger &a, const BigInteger &b);
+	/* See the comment on BigUnsigned::divideWithRemainder.  Semantics
+	 * differ from those of primitive integers when negatives and/or zeros
+	 * are involved. */
+	void divideWithRemainder(const BigInteger &b, BigInteger &q);
+	void negate(const BigInteger &a);
+	
+	/* Bitwise operators are not provided for BigIntegers.  Use
+	 * getMagnitude to get the magnitude and operate on that instead. */
+
+	BigInteger operator +(const BigInteger &x) const;
+	BigInteger operator -(const BigInteger &x) const;
+	BigInteger operator *(const BigInteger &x) const;
+	BigInteger operator /(const BigInteger &x) const;
+	BigInteger operator %(const BigInteger &x) const;
+	BigInteger operator -() const;
+
+	void operator +=(const BigInteger &x);
+	void operator -=(const BigInteger &x);
+	void operator *=(const BigInteger &x);
+	void operator /=(const BigInteger &x);
+	void operator %=(const BigInteger &x);
+	void flipSign();
+
+	// INCREMENT/DECREMENT OPERATORS
+	void operator ++(   );
+	void operator ++(int);
+	void operator --(   );
+	void operator --(int);
+};
+
+// NORMAL OPERATORS
+/* These create an object to hold the result and invoke
+ * the appropriate put-here operation on it, passing
+ * this and x.  The new object is then returned. */
+inline BigInteger BigInteger::operator +(const BigInteger &x) const {
+	BigInteger ans;
+	ans.add(*this, x);
+	return ans;
+}
+inline BigInteger BigInteger::operator -(const BigInteger &x) const {
+	BigInteger ans;
+	ans.subtract(*this, x);
+	return ans;
+}
+inline BigInteger BigInteger::operator *(const BigInteger &x) const {
+	BigInteger ans;
+	ans.multiply(*this, x);
+	return ans;
+}
+inline BigInteger BigInteger::operator /(const BigInteger &x) const {
+	if (x.isZero()) throw "BigInteger::operator /: division by zero";
+	BigInteger q, r;
+	r = *this;
+	r.divideWithRemainder(x, q);
+	return q;
+}
+inline BigInteger BigInteger::operator %(const BigInteger &x) const {
+	if (x.isZero()) throw "BigInteger::operator %: division by zero";
+	BigInteger q, r;
+	r = *this;
+	r.divideWithRemainder(x, q);
+	return r;
+}
+inline BigInteger BigInteger::operator -() const {
+	BigInteger ans;
+	ans.negate(*this);
+	return ans;
+}
+
+/*
+ * ASSIGNMENT OPERATORS
+ * 
+ * Now the responsibility for making a temporary copy if necessary
+ * belongs to the put-here operations.  See Assignment Operators in
+ * BigUnsigned.hh.
+ */
+inline void BigInteger::operator +=(const BigInteger &x) {
+	add(*this, x);
+}
+inline void BigInteger::operator -=(const BigInteger &x) {
+	subtract(*this, x);
+}
+inline void BigInteger::operator *=(const BigInteger &x) {
+	multiply(*this, x);
+}
+inline void BigInteger::operator /=(const BigInteger &x) {
+	if (x.isZero()) throw "BigInteger::operator /=: division by zero";
+	/* The following technique is slightly faster than copying *this first
+	 * when x is large. */
+	BigInteger q;
+	divideWithRemainder(x, q);
+	// *this contains the remainder, but we overwrite it with the quotient.
+	*this = q;
+}
+inline void BigInteger::operator %=(const BigInteger &x) {
+	if (x.isZero()) throw "BigInteger::operator %=: division by zero";
+	BigInteger q;
+	// Mods *this by x.  Don't care about quotient left in q.
+	divideWithRemainder(x, q);
+}
+// This one is trivial
+inline void BigInteger::flipSign() {
+	sign = Sign(-sign);
+}
+
+#endif
diff --git a/bigint/BigIntegerAlgorithms.cc b/bigint/BigIntegerAlgorithms.cc
new file mode 100644
index 0000000..7edebda
--- /dev/null
+++ b/bigint/BigIntegerAlgorithms.cc
@@ -0,0 +1,70 @@
+#include "BigIntegerAlgorithms.hh"
+
+BigUnsigned gcd(BigUnsigned a, BigUnsigned b) {
+	BigUnsigned trash;
+	// Neat in-place alternating technique.
+	for (;;) {
+		if (b.isZero())
+			return a;
+		a.divideWithRemainder(b, trash);
+		if (a.isZero())
+			return b;
+		b.divideWithRemainder(a, trash);
+	}
+}
+
+void extendedEuclidean(BigInteger m, BigInteger n,
+		BigInteger &g, BigInteger &r, BigInteger &s) {
+	if (&g == &r || &g == &s || &r == &s)
+		throw "BigInteger extendedEuclidean: Outputs are aliased";
+	BigInteger r1(1), s1(0), r2(0), s2(1), q;
+	/* Invariants:
+	 * r1*m(orig) + s1*n(orig) == m(current)
+	 * r2*m(orig) + s2*n(orig) == n(current) */
+	for (;;) {
+		if (n.isZero()) {
+			r = r1; s = s1; g = m;
+			return;
+		}
+		// Subtract q times the second invariant from the first invariant.
+		m.divideWithRemainder(n, q);
+		r1 -= q*r2; s1 -= q*s2;
+
+		if (m.isZero()) {
+			r = r2; s = s2; g = n;
+			return;
+		}
+		// Subtract q times the first invariant from the second invariant.
+		n.divideWithRemainder(m, q);
+		r2 -= q*r1; s2 -= q*s1;
+	}
+}
+
+BigUnsigned modinv(const BigInteger &x, const BigUnsigned &n) {
+	BigInteger g, r, s;
+	extendedEuclidean(x, n, g, r, s);
+	if (g == 1)
+		// r*x + s*n == 1, so r*x === 1 (mod n), so r is the answer.
+		return (r % n).getMagnitude(); // (r % n) will be nonnegative
+	else
+		throw "BigInteger modinv: x and n have a common factor";
+}
+
+BigUnsigned modexp(const BigInteger &base, const BigUnsigned &exponent,
+		const BigUnsigned &modulus) {
+	BigUnsigned ans = 1, base2 = (base % modulus).getMagnitude();
+	BigUnsigned::Index i = exponent.bitLength();
+	// For each bit of the exponent, most to least significant...
+	while (i > 0) {
+		i--;
+		// Square.
+		ans *= ans;
+		ans %= modulus;
+		// And multiply if the bit is a 1.
+		if (exponent.getBit(i)) {
+			ans *= base2;
+			ans %= modulus;
+		}
+	}
+	return ans;
+}
diff --git a/bigint/BigIntegerAlgorithms.hh b/bigint/BigIntegerAlgorithms.hh
new file mode 100644
index 0000000..b1dd943
--- /dev/null
+++ b/bigint/BigIntegerAlgorithms.hh
@@ -0,0 +1,25 @@
+#ifndef BIGINTEGERALGORITHMS_H
+#define BIGINTEGERALGORITHMS_H
+
+#include "BigInteger.hh"
+
+/* Some mathematical algorithms for big integers.
+ * This code is new and, as such, experimental. */
+
+// Returns the greatest common divisor of a and b.
+BigUnsigned gcd(BigUnsigned a, BigUnsigned b);
+
+/* Extended Euclidean algorithm.
+ * Given m and n, finds gcd g and numbers r, s such that r*m + s*n == g. */
+void extendedEuclidean(BigInteger m, BigInteger n,
+		BigInteger &g, BigInteger &r, BigInteger &s);
+
+/* Returns the multiplicative inverse of x modulo n, or throws an exception if
+ * they have a common factor. */
+BigUnsigned modinv(const BigInteger &x, const BigUnsigned &n);
+
+// Returns (base ^ exponent) % modulus.
+BigUnsigned modexp(const BigInteger &base, const BigUnsigned &exponent,
+		const BigUnsigned &modulus);
+
+#endif
diff --git a/bigint/BigIntegerLibrary.hh b/bigint/BigIntegerLibrary.hh
new file mode 100644
index 0000000..2a0ebee
--- /dev/null
+++ b/bigint/BigIntegerLibrary.hh
@@ -0,0 +1,8 @@
+// This header file includes all of the library header files.
+
+#include "NumberlikeArray.hh"
+#include "BigUnsigned.hh"
+#include "BigInteger.hh"
+#include "BigIntegerAlgorithms.hh"
+#include "BigUnsignedInABase.hh"
+#include "BigIntegerUtils.hh"
diff --git a/bigint/BigIntegerUtils.cc b/bigint/BigIntegerUtils.cc
new file mode 100644
index 0000000..44073af
--- /dev/null
+++ b/bigint/BigIntegerUtils.cc
@@ -0,0 +1,50 @@
+#include "BigIntegerUtils.hh"
+#include "BigUnsignedInABase.hh"
+
+std::string bigUnsignedToString(const BigUnsigned &x) {
+	return std::string(BigUnsignedInABase(x, 10));
+}
+
+std::string bigIntegerToString(const BigInteger &x) {
+	return (x.getSign() == BigInteger::negative)
+		? (std::string("-") + bigUnsignedToString(x.getMagnitude()))
+		: (bigUnsignedToString(x.getMagnitude()));
+}
+
+BigUnsigned stringToBigUnsigned(const std::string &s) {
+	return BigUnsigned(BigUnsignedInABase(s, 10));
+}
+
+BigInteger stringToBigInteger(const std::string &s) {
+	// Recognize a sign followed by a BigUnsigned.
+	return (s[0] == '-') ? BigInteger(stringToBigUnsigned(s.substr(1, s.length() - 1)), BigInteger::negative)
+		: (s[0] == '+') ? BigInteger(stringToBigUnsigned(s.substr(1, s.length() - 1)))
+		: BigInteger(stringToBigUnsigned(s));
+}
+
+std::ostream &operator <<(std::ostream &os, const BigUnsigned &x) {
+	BigUnsignedInABase::Base base;
+	long osFlags = os.flags();
+	if (osFlags & os.dec)
+		base = 10;
+	else if (osFlags & os.hex) {
+		base = 16;
+		if (osFlags & os.showbase)
+			os << "0x";
+	} else if (osFlags & os.oct) {
+		base = 8;
+		if (osFlags & os.showbase)
+			os << '0';
+	} else
+		throw "std::ostream << BigUnsigned: Could not determine the desired base from output-stream flags";
+	std::string s = std::string(BigUnsignedInABase(x, base));
+	os << s;
+	return os;
+}
+
+std::ostream &operator <<(std::ostream &os, const BigInteger &x) {
+	if (x.getSign() == BigInteger::negative)
+		os << '-';
+	os << x.getMagnitude();
+	return os;
+}
diff --git a/bigint/BigIntegerUtils.hh b/bigint/BigIntegerUtils.hh
new file mode 100644
index 0000000..c815b5d
--- /dev/null
+++ b/bigint/BigIntegerUtils.hh
@@ -0,0 +1,72 @@
+#ifndef BIGINTEGERUTILS_H
+#define BIGINTEGERUTILS_H
+
+#include "BigInteger.hh"
+#include <string>
+#include <iostream>
+
+/* This file provides:
+ * - Convenient std::string <-> BigUnsigned/BigInteger conversion routines
+ * - std::ostream << operators for BigUnsigned/BigInteger */
+
+// std::string conversion routines.  Base 10 only.
+std::string bigUnsignedToString(const BigUnsigned &x);
+std::string bigIntegerToString(const BigInteger &x);
+BigUnsigned stringToBigUnsigned(const std::string &s);
+BigInteger stringToBigInteger(const std::string &s);
+
+// Creates a BigInteger from data such as `char's; read below for details.
+template <class T>
+BigInteger dataToBigInteger(const T* data, BigInteger::Index length, BigInteger::Sign sign);
+
+// Outputs x to os, obeying the flags `dec', `hex', `bin', and `showbase'.
+std::ostream &operator <<(std::ostream &os, const BigUnsigned &x);
+
+// Outputs x to os, obeying the flags `dec', `hex', `bin', and `showbase'.
+// My somewhat arbitrary policy: a negative sign comes before a base indicator (like -0xFF).
+std::ostream &operator <<(std::ostream &os, const BigInteger &x);
+
+// BEGIN TEMPLATE DEFINITIONS.
+
+/*
+ * Converts binary data to a BigInteger.
+ * Pass an array `data', its length, and the desired sign.
+ *
+ * Elements of `data' may be of any type `T' that has the following
+ * two properties (this includes almost all integral types):
+ *
+ * (1) `sizeof(T)' correctly gives the amount of binary data in one
+ * value of `T' and is a factor of `sizeof(Blk)'.
+ *
+ * (2) When a value of `T' is casted to a `Blk', the low bytes of
+ * the result contain the desired binary data.
+ */
+template <class T>
+BigInteger dataToBigInteger(const T* data, BigInteger::Index length, BigInteger::Sign sign) {
+	// really ceiling(numBytes / sizeof(BigInteger::Blk))
+	unsigned int pieceSizeInBits = 8 * sizeof(T);
+	unsigned int piecesPerBlock = sizeof(BigInteger::Blk) / sizeof(T);
+	unsigned int numBlocks = (length + piecesPerBlock - 1) / piecesPerBlock;
+
+	// Allocate our block array
+	BigInteger::Blk *blocks = new BigInteger::Blk[numBlocks];
+
+	BigInteger::Index blockNum, pieceNum, pieceNumHere;
+
+	// Convert
+	for (blockNum = 0, pieceNum = 0; blockNum < numBlocks; blockNum++) {
+		BigInteger::Blk curBlock = 0;
+		for (pieceNumHere = 0; pieceNumHere < piecesPerBlock && pieceNum < length;
+			pieceNumHere++, pieceNum++)
+			curBlock |= (BigInteger::Blk(data[pieceNum]) << (pieceSizeInBits * pieceNumHere));
+		blocks[blockNum] = curBlock;
+	}
+
+	// Create the BigInteger.
+	BigInteger x(blocks, numBlocks, sign);
+
+	delete [] blocks;
+	return x;
+}
+
+#endif
diff --git a/bigint/BigUnsigned.cc b/bigint/BigUnsigned.cc
new file mode 100644
index 0000000..d7f9889
--- /dev/null
+++ b/bigint/BigUnsigned.cc
@@ -0,0 +1,697 @@
+#include "BigUnsigned.hh"
+
+// Memory management definitions have moved to the bottom of NumberlikeArray.hh.
+
+// The templates used by these constructors and converters are at the bottom of
+// BigUnsigned.hh.
+
+BigUnsigned::BigUnsigned(unsigned long  x) { initFromPrimitive      (x); }
+BigUnsigned::BigUnsigned(unsigned int   x) { initFromPrimitive      (x); }
+BigUnsigned::BigUnsigned(unsigned short x) { initFromPrimitive      (x); }
+BigUnsigned::BigUnsigned(         long  x) { initFromSignedPrimitive(x); }
+BigUnsigned::BigUnsigned(         int   x) { initFromSignedPrimitive(x); }
+BigUnsigned::BigUnsigned(         short x) { initFromSignedPrimitive(x); }
+
+unsigned long  BigUnsigned::toUnsignedLong () const { return convertToPrimitive      <unsigned long >(); }
+unsigned int   BigUnsigned::toUnsignedInt  () const { return convertToPrimitive      <unsigned int  >(); }
+unsigned short BigUnsigned::toUnsignedShort() const { return convertToPrimitive      <unsigned short>(); }
+long           BigUnsigned::toLong         () const { return convertToSignedPrimitive<         long >(); }
+int            BigUnsigned::toInt          () const { return convertToSignedPrimitive<         int  >(); }
+short          BigUnsigned::toShort        () const { return convertToSignedPrimitive<         short>(); }
+
+// BIT/BLOCK ACCESSORS
+
+void BigUnsigned::setBlock(Index i, Blk newBlock) {
+	if (newBlock == 0) {
+		if (i < len) {
+			blk[i] = 0;
+			zapLeadingZeros();
+		}
+		// If i >= len, no effect.
+	} else {
+		if (i >= len) {
+			// The nonzero block extends the number.
+			allocateAndCopy(i+1);
+			// Zero any added blocks that we aren't setting.
+			for (Index j = len; j < i; j++)
+				blk[j] = 0;
+			len = i+1;
+		}
+		blk[i] = newBlock;
+	}
+}
+
+/* Evidently the compiler wants BigUnsigned:: on the return type because, at
+ * that point, it hasn't yet parsed the BigUnsigned:: on the name to get the
+ * proper scope. */
+BigUnsigned::Index BigUnsigned::bitLength() const {
+	if (isZero())
+		return 0;
+	else {
+		Blk leftmostBlock = getBlock(len - 1);