bignum: deprecate by moving into exp directory

R=rsc
CC=golang-dev
http://codereview.appspot.com/1211047

16 files changed, 2802 insertions, 8 deletions

diff --git a/src/pkg/exp/bignum/Makefile b/src/pkg/exp/bignum/Makefile
new file mode 100644
index 000000000..064cf1eb9
--- /dev/null
+++ b/src/pkg/exp/bignum/Makefile
@@ -0,0 +1,14 @@
+# Copyright 2009 The Go Authors. All rights reserved.
+# Use of this source code is governed by a BSD-style
+# license that can be found in the LICENSE file.
+
+include ../../../Make.$(GOARCH)
+
+TARG=exp/bignum
+GOFILES=\
+	arith.go\
+	bignum.go\
+	integer.go\
+	rational.go\
+
+include ../../../Make.pkg
diff --git a/src/pkg/exp/bignum/arith.go b/src/pkg/exp/bignum/arith.go
new file mode 100644
index 000000000..aa65dbd7a
--- /dev/null
+++ b/src/pkg/exp/bignum/arith.go
@@ -0,0 +1,121 @@
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Fast versions of the routines in this file are in fast.arith.s.
+// Simply replace this file with arith.s (renamed from fast.arith.s)
+// and the bignum package will build and run on a platform that
+// supports the assembly routines.
+
+package bignum
+
+import "unsafe"
+
+// z1<<64 + z0 = x*y
+func Mul128(x, y uint64) (z1, z0 uint64) {
+	// Split x and y into 2 halfwords each, multiply
+	// the halfwords separately while avoiding overflow,
+	// and return the product as 2 words.
+
+	const (
+		W  = uint(unsafe.Sizeof(x)) * 8
+		W2 = W / 2
+		B2 = 1 << W2
+		M2 = B2 - 1
+	)
+
+	if x < y {
+		x, y = y, x
+	}
+
+	if x < B2 {
+		// y < B2 because y <= x
+		// sub-digits of x and y are (0, x) and (0, y)
+		// z = z[0] = x*y
+		z0 = x * y
+		return
+	}
+
+	if y < B2 {
+		// sub-digits of x and y are (x1, x0) and (0, y)
+		// x = (x1*B2 + x0)
+		// y = (y1*B2 + y0)
+		x1, x0 := x>>W2, x&M2
+
+		// x*y = t2*B2*B2 + t1*B2 + t0
+		t0 := x0 * y
+		t1 := x1 * y
+
+		// compute result digits but avoid overflow
+		// z = z[1]*B + z[0] = x*y
+		z0 = t1<<W2 + t0
+		z1 = (t1 + t0>>W2) >> W2
+		return
+	}
+
+	// general case
+	// sub-digits of x and y are (x1, x0) and (y1, y0)
+	// x = (x1*B2 + x0)
+	// y = (y1*B2 + y0)
+	x1, x0 := x>>W2, x&M2
+	y1, y0 := y>>W2, y&M2
+
+	// x*y = t2*B2*B2 + t1*B2 + t0
+	t0 := x0 * y0
+	t1 := x1*y0 + x0*y1
+	t2 := x1 * y1
+
+	// compute result digits but avoid overflow
+	// z = z[1]*B + z[0] = x*y
+	z0 = t1<<W2 + t0
+	z1 = t2 + (t1+t0>>W2)>>W2
+	return
+}
+
+
+// z1<<64 + z0 = x*y + c
+func MulAdd128(x, y, c uint64) (z1, z0 uint64) {
+	// Split x and y into 2 halfwords each, multiply
+	// the halfwords separately while avoiding overflow,
+	// and return the product as 2 words.
+
+	const (
+		W  = uint(unsafe.Sizeof(x)) * 8
+		W2 = W / 2
+		B2 = 1 << W2
+		M2 = B2 - 1
+	)
+
+	// TODO(gri) Should implement special cases for faster execution.
+
+	// general case
+	// sub-digits of x, y, and c are (x1, x0), (y1, y0), (c1, c0)
+	// x = (x1*B2 + x0)
+	// y = (y1*B2 + y0)
+	x1, x0 := x>>W2, x&M2
+	y1, y0 := y>>W2, y&M2
+	c1, c0 := c>>W2, c&M2
+
+	// x*y + c = t2*B2*B2 + t1*B2 + t0
+	t0 := x0*y0 + c0
+	t1 := x1*y0 + x0*y1 + c1
+	t2 := x1 * y1
+
+	// compute result digits but avoid overflow
+	// z = z[1]*B + z[0] = x*y
+	z0 = t1<<W2 + t0
+	z1 = t2 + (t1+t0>>W2)>>W2
+	return
+}
+
+
+// q = (x1<<64 + x0)/y + r
+func Div128(x1, x0, y uint64) (q, r uint64) {
+	if x1 == 0 {
+		q, r = x0/y, x0%y
+		return
+	}
+
+	// TODO(gri) Implement general case.
+	panic("Div128 not implemented for x > 1<<64-1")
+}
diff --git a/src/pkg/exp/bignum/arith_amd64.s b/src/pkg/exp/bignum/arith_amd64.s
new file mode 100644
index 000000000..37d5a30de
--- /dev/null
+++ b/src/pkg/exp/bignum/arith_amd64.s
@@ -0,0 +1,41 @@
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// This file provides fast assembly versions
+// of the routines in arith.go.
+
+// func Mul128(x, y uint64) (z1, z0 uint64)
+// z1<<64 + z0 = x*y
+//
+TEXT ·Mul128(SB),7,$0
+	MOVQ a+0(FP), AX
+	MULQ a+8(FP)
+	MOVQ DX, a+16(FP)
+	MOVQ AX, a+24(FP)
+	RET
+
+
+// func MulAdd128(x, y, c uint64) (z1, z0 uint64)
+// z1<<64 + z0 = x*y + c
+//
+TEXT ·MulAdd128(SB),7,$0
+	MOVQ a+0(FP), AX
+	MULQ a+8(FP)
+	ADDQ a+16(FP), AX
+	ADCQ $0, DX
+	MOVQ DX, a+24(FP)
+	MOVQ AX, a+32(FP)
+	RET
+
+
+// func Div128(x1, x0, y uint64) (q, r uint64)
+// q = (x1<<64 + x0)/y + r
+//
+TEXT ·Div128(SB),7,$0
+	MOVQ a+0(FP), DX
+	MOVQ a+8(FP), AX
+	DIVQ a+16(FP)
+	MOVQ AX, a+24(FP)
+	MOVQ DX, a+32(FP)
+	RET
diff --git a/src/pkg/exp/bignum/bignum.go b/src/pkg/exp/bignum/bignum.go
new file mode 100644
index 000000000..485583199
--- /dev/null
+++ b/src/pkg/exp/bignum/bignum.go
@@ -0,0 +1,1024 @@
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// A package for arbitrary precision arithmethic.
+// It implements the following numeric types:
+//
+//	- Natural	unsigned integers
+//	- Integer	signed integers
+//	- Rational	rational numbers
+//
+// This package has been designed for ease of use but the functions it provides
+// are likely to be quite slow. It may be deprecated eventually. Use package
+// big instead, if possible.
+//
+package bignum
+
+import (
+	"fmt"
+)
+
+// TODO(gri) Complete the set of in-place operations.
+
+// ----------------------------------------------------------------------------
+// Internal representation
+//
+// A natural number of the form
+//
+//   x = x[n-1]*B^(n-1) + x[n-2]*B^(n-2) + ... + x[1]*B + x[0]
+//
+// with 0 <= x[i] < B and 0 <= i < n is stored in a slice of length n,
+// with the digits x[i] as the slice elements.
+//
+// A natural number is normalized if the slice contains no leading 0 digits.
+// During arithmetic operations, denormalized values may occur but are
+// always normalized before returning the final result. The normalized
+// representation of 0 is the empty slice (length = 0).
+//
+// The operations for all other numeric types are implemented on top of
+// the operations for natural numbers.
+//
+// The base B is chosen as large as possible on a given platform but there
+// are a few constraints besides the size of the largest unsigned integer
+// type available:
+//
+// 1) To improve conversion speed between strings and numbers, the base B
+//    is chosen such that division and multiplication by 10 (for decimal
+//    string representation) can be done without using extended-precision
+//    arithmetic. This makes addition, subtraction, and conversion routines
+//    twice as fast. It requires a ``buffer'' of 4 bits per operand digit.
+//    That is, the size of B must be 4 bits smaller then the size of the
+//    type (digit) in which these operations are performed. Having this
+//    buffer also allows for trivial (single-bit) carry computation in
+//    addition and subtraction (optimization suggested by Ken Thompson).
+//
+// 2) Long division requires extended-precision (2-digit) division per digit.
+//    Instead of sacrificing the largest base type for all other operations,
+//    for division the operands are unpacked into ``half-digits'', and the
+//    results are packed again. For faster unpacking/packing, the base size
+//    in bits must be even.
+
+type (
+	digit  uint64
+	digit2 uint32 // half-digits for division
+)
+
+
+const (
+	logW = 64          // word width
+	logH = 4           // bits for a hex digit (= small number)
+	logB = logW - logH // largest bit-width available
+
+	// half-digits
+	_W2 = logB / 2 // width
+	_B2 = 1 << _W2 // base
+	_M2 = _B2 - 1  // mask
+
+	// full digits
+	_W = _W2 * 2 // width
+	_B = 1 << _W // base
+	_M = _B - 1  // mask
+)
+
+
+// ----------------------------------------------------------------------------
+// Support functions
+
+func assert(p bool) {
+	if !p {
+		panic("assert failed")
+	}
+}
+
+
+func isSmall(x digit) bool { return x < 1<<logH }
+
+
+// For debugging. Keep around.
+/*
+func dump(x Natural) {
+	print("[", len(x), "]");
+	for i := len(x) - 1; i >= 0; i-- {
+		print(" ", x[i]);
+	}
+	println();
+}
+*/
+
+
+// ----------------------------------------------------------------------------
+// Natural numbers
+
+// Natural represents an unsigned integer value of arbitrary precision.
+//
+type Natural []digit
+
+
+// Nat creates a small natural number with value x.
+//
+func Nat(x uint64) Natural {
+	if x == 0 {
+		return nil // len == 0
+	}
+
+	// single-digit values
+	// (note: cannot re-use preallocated values because
+	//        the in-place operations may overwrite them)
+	if x < _B {
+		return Natural{digit(x)}
+	}
+
+	// compute number of digits required to represent x
+	// (this is usually 1 or 2, but the algorithm works
+	// for any base)
+	n := 0
+	for t := x; t > 0; t >>= _W {
+		n++
+	}
+
+	// split x into digits
+	z := make(Natural, n)
+	for i := 0; i < n; i++ {
+		z[i] = digit(x & _M)
+		x >>= _W
+	}
+
+	return z
+}
+
+
+// Value returns the lowest 64bits of x.
+//
+func (x Natural) Value() uint64 {
+	// single-digit values
+	n := len(x)
+	switch n {
+	case 0:
+		return 0
+	case 1:
+		return uint64(x[0])
+	}
+
+	// multi-digit values
+	// (this is usually 1 or 2, but the algorithm works
+	// for any base)
+	z := uint64(0)
+	s := uint(0)
+	for i := 0; i < n && s < 64; i++ {
+		z += uint64(x[i]) << s
+		s += _W
+	}
+
+	return z
+}
+
+
+// Predicates
+
+// IsEven returns true iff x is divisible by 2.
+//
+func (x Natural) IsEven() bool { return len(x) == 0 || x[0]&1 == 0 }
+
+
+// IsOdd returns true iff x is not divisible by 2.
+//
+func (x Natural) IsOdd() bool { return len(x) > 0 && x[0]&1 != 0 }
+
+
+// IsZero returns true iff x == 0.
+//
+func (x Natural) IsZero() bool { return len(x) == 0 }
+
+
+// Operations
+//
+// Naming conventions
+//
+// c      carry
+// x, y   operands
+// z      result
+// n, m   len(x), len(y)
+
+func normalize(x Natural) Natural {
+	n := len(x)
+	for n > 0 && x[n-1] == 0 {
+		n--
+	}
+	return x[0:n]
+}
+
+
+// nalloc returns a Natural of n digits. If z is large
+// enough, z is resized and returned. Otherwise, a new
+// Natural is allocated.
+//
+func nalloc(z Natural, n int) Natural {
+	size := n
+	if size <= 0 {
+		size = 4
+	}
+	if size <= cap(z) {
+		return z[0:n]
+	}
+	return make(Natural, n, size)
+}
+
+
+// Nadd sets *zp to the sum x + y.
+// *zp may be x or y.
+//
+func Nadd(zp *Natural, x, y Natural) {
+	n := len(x)
+	m := len(y)
+	if n < m {
+		Nadd(zp, y, x)
+		return
+	}
+
+	z := nalloc(*zp, n+1)
+	c := digit(0)
+	i := 0
+	for i < m {
+		t := c + x[i] + y[i]
+		c, z[i] = t>>_W, t&_M
+		i++
+	}
+	for i < n {
+		t := c + x[i]
+		c, z[i] = t>>_W, t&_M
+		i++
+	}
+	if c != 0 {
+		z[i] = c
+		i++
+	}
+	*zp = z[0:i]
+}
+
+
+// Add returns the sum z = x + y.
+//
+func (x Natural) Add(y Natural) Natural {
+	var z Natural
+	Nadd(&z, x, y)
+	return z
+}
+
+
+// Nsub sets *zp to the difference x - y for x >= y.
+// If x < y, an underflow run-time error occurs (use Cmp to test if x >= y).
+// *zp may be x or y.
+//
+func Nsub(zp *Natural, x, y Natural) {
+	n := len(x)
+	m := len(y)
+	if n < m {
+		panic("underflow")
+	}
+
+	z := nalloc(*zp, n)
+	c := digit(0)
+	i := 0
+	for i < m {
+		t := c + x[i] - y[i]
+		c, z[i] = digit(int64(t)>>_W), t&_M // requires arithmetic shift!
+		i++
+	}
+	for i < n {
+		t := c + x[i]
+		c, z[i] = digit(int64(t)>>_W), t&_M // requires arithmetic shift!
+		i++
+	}
+	if int64(c) < 0 {
+		panic("underflow")
+	}
+	*zp = normalize(z)
+}
+
+
+// Sub returns the difference x - y for x >= y.
+// If x < y, an underflow run-time error occurs (use Cmp to test if x >= y).
+//
+func (x Natural) Sub(y Natural) Natural {
+	var z Natural
+	Nsub(&z, x, y)
+	return z
+}
+
+
+// Returns z1 = (x*y + c) div B, z0 = (x*y + c) mod B.
+//
+func muladd11(x, y, c digit) (digit, digit) {
+	z1, z0 := MulAdd128(uint64(x), uint64(y), uint64(c))
+	return digit(z1<<(64-logB) | z0>>logB), digit(z0 & _M)
+}
+
+
+func mul1(z, x Natural, y digit) (c digit) {
+	for i := 0; i < len(x); i++ {
+		c, z[i] = muladd11(x[i], y, c)
+	}
+	return
+}
+
+
+// Nscale sets *z to the scaled value (*z) * d.
+//
+func Nscale(z *Natural, d uint64) {
+	switch {
+	case d == 0:
+		*z = Nat(0)
+		return
+	case d == 1:
+		return
+	case d >= _B:
+		*z = z.Mul1(d)
+		return
+	}
+
+	c := mul1(*z, *z, digit(d))
+
+	if c != 0 {
+		n := len(*z)
+		if n >= cap(*z) {
+			zz := make(Natural, n+1)
+			for i, d := range *z {
+				zz[i] = d
+			}
+			*z = zz
+		} else {
+			*z = (*z)[0 : n+1]
+		}
+		(*z)[n] = c
+	}
+}
+
+
+// Computes x = x*d + c for small d's.
+//
+func muladd1(x Natural, d, c digit) Natural {
+	assert(isSmall(d-1) && isSmall(c))
+	n := len(x)
+	z := make(Natural, n+1)
+
+	for i := 0; i < n; i++ {
+		t := c + x[i]*d
+		c, z[i] = t>>_W, t&_M
+	}
+	z[n] = c
+
+	return normalize(z)
+}
+
+
+// Mul1 returns the product x * d.
+//
+func (x Natural) Mul1(d uint64) Natural {
+	switch {
+	case d == 0:
+		return Nat(0)
+	case d == 1:
+		return x
+	case isSmall(digit(d)):
+		muladd1(x, digit(d), 0)
+	case d >= _B:
+		return x.Mul(Nat(d))
+	}
+
+	z := make(Natural, len(x)+1)
+	c := mul1(z, x, digit(d))
+	z[len(x)] = c
+	return normalize(z)
+}
+
+
+// Mul returns the product x * y.
+//
+func (x Natural) Mul(y Natural) Natural {
+	n := len(x)
+	m := len(y)
+	if n < m {
+		return y.Mul(x)
+	}
+
+	if m == 0 {
+		return Nat(0)
+	}
+
+	if m == 1 && y[0] < _B {
+		return x.Mul1(uint64(y[0]))
+	}
+
+	z := make(Natural, n+m)
+	for j := 0; j < m; j++ {
+		d := y[j]
+		if d != 0 {
+			c := digit(0)
+			for i := 0; i < n; i++ {
+				c, z[i+j] = muladd11(x[i], d, z[i+j]+c)
+			}
+			z[n+j] = c
+		}
+	}
+
+	return normalize(z)
+}
+
+
+// DivMod needs multi-precision division, which is not available if digit
+// is already using the largest uint size. Instead, unpack each operand
+// into operands with twice as many digits of half the size (digit2), do
+// DivMod, and then pack the results again.
+
+func unpack(x Natural) []digit2 {
+	n := len(x)
+	z := make([]digit2, n*2+1) // add space for extra digit (used by DivMod)
+	for i := 0; i < n; i++ {
+		t := x[i]
+		z[i*2] = digit2(t & _M2)
+		z[i*2+1] = digit2(t >> _W2 & _M2)
+	}
+
+	// normalize result
+	k := 2 * n
+	for k > 0 && z[k-1] == 0 {
+		k--
+	}
+	return z[0:k] // trim leading 0's
+}
+
+
+func pack(x []digit2) Natural {
+	n := (len(x) + 1) / 2
+	z := make(Natural, n)
+	if len(x)&1 == 1 {
+		// handle odd len(x)
+		n--
+		z[n] = digit(x[n*2])
+	}
+	for i := 0; i < n; i++ {
+		z[i] = digit(x[i*2+1])<<_W2 | digit(x[i*2])
+	}
+	return normalize(z)
+}
+
+
+func mul21(z, x []digit2, y digit2) digit2 {
+	c := digit(0)
+	f := digit(y)
+	for i := 0; i < len(x); i++ {
+		t := c + digit(x[i])*f
+		c, z[i] = t>>_W2, digit2(t&_M2)
+	}
+	return digit2(c)
+}
+
+
+func div21(z, x []digit2, y digit2) digit2 {
+	c := digit(0)
+	d := digit(y)
+	for i := len(x) - 1; i >= 0; i-- {
+		t := c<<_W2 + digit(x[i])
+		c, z[i] = t%d, digit2(t/d)
+	}
+	return digit2(c)
+}
+
+
+// divmod returns q and r with x = y*q + r and 0 <= r < y.
+// x and y are destroyed in the process.
+//
+// The algorithm used here is based on 1). 2) describes the same algorithm
+// in C. A discussion and summary of the relevant theorems can be found in
+// 3). 3) also describes an easier way to obtain the trial digit - however
+// it relies on tripple-precision arithmetic which is why Knuth's method is
+// used here.
+//
+// 1) D. Knuth, The Art of Computer Programming. Volume 2. Seminumerical
+//    Algorithms. Addison-Wesley, Reading, 1969.
+//    (Algorithm D, Sec. 4.3.1)
+//
+// 2) Henry S. Warren, Jr., Hacker's Delight. Addison-Wesley, 2003.
+//    (9-2 Multiword Division, p.140ff)
+//
+// 3) P. Brinch Hansen, ``Multiple-length division revisited: A tour of the
+//    minefield''. Software - Practice and Experience 24, (June 1994),
+//    579-601. John Wiley & Sons, Ltd.
+
+func divmod(x, y []digit2) ([]digit2, []digit2) {
+	n := len(x)
+	m := len(y)
+	if m == 0 {
+		panic("division by zero")
+	}
+	assert(n+1 <= cap(x)) // space for one extra digit
+	x = x[0 : n+1]
+	assert(x[n] == 0)
+
+	if m == 1 {
+		// division by single digit
+		// result is shifted left by 1 in place!
+		x[0] = div21(x[1:n+1], x[0:n], y[0])
+
+	} else if m > n {
+		// y > x => quotient = 0, remainder = x
+		// TODO in this case we shouldn't even unpack x and y
+		m = n
+
+	} else {
+		// general case
+		assert(2 <= m && m <= n)
+
+		// normalize x and y
+		// TODO Instead of multiplying, it would be sufficient to
+		//      shift y such that the normalization condition is
+		//      satisfied (as done in Hacker's Delight).
+		f := _B2 / (digit(y[m-1]) + 1)
+		if f != 1 {
+			mul21(x, x, digit2(f))
+			mul21(y, y, digit2(f))
+		}
+		assert(_B2/2 <= y[m-1] && y[m-1] < _B2) // incorrect scaling
+
+		y1, y2 := digit(y[m-1]), digit(y[m-2])
+		for i := n - m; i >= 0; i-- {
+			k := i + m
+
+			// compute trial digit (Knuth)
+			var q digit
+			{
+				x0, x1, x2 := digit(x[k]), digit(x[k-1]), digit(x[k-2])
+				if x0 != y1 {
+					q = (x0<<_W2 + x1) / y1
+				} else {
+					q = _B2 - 1
+				}
+				for y2*q > (x0<<_W2+x1-y1*q)<<_W2+x2 {
+					q--
+				}
+			}
+
+			// subtract y*q
+			c := digit(0)
+			for j := 0; j < m; j++ {
+				t := c + digit(x[i+j]) - digit(y[j])*q
+				c, x[i+j] = digit(int64(t)>>_W2), digit2(t&_M2) // requires arithmetic shift!
+			}
+			x[k] = digit2((c + digit(x[k])) & _M2)
+
+			// correct if trial digit was too large
+			if x[k] != 0 {
+				// add y
+				c := digit(0)
+				for j := 0; j < m; j++ {
+					t := c + digit(x[i+j]) + digit(y[j])
+					c, x[i+j] = t>>_W2, digit2(t&_M2)
+				}
+				x[k] = digit2((c + digit(x[k])) & _M2)
+				assert(x[k] == 0)
+				// correct trial digit
+				q--
+			}
+
+			x[k] = digit2(q)
+		}
+
+		// undo normalization for remainder
+		if f != 1 {
+			c := div21(x[0:m], x[0:m], digit2(f))
+			assert(c == 0)
+		}
+	}
+
+	return x[m : n+1], x[0:m]
+}
+
+
+// Div returns the quotient q = x / y for y > 0,
+// with x = y*q + r and 0 <= r < y.
+// If y == 0, a division-by-zero run-time error occurs.
+//
+func (x Natural) Div(y Natural) Natural {
+	q, _ := divmod(unpack(x), unpack(y))
+	return pack(q)
+}
+
+
+// Mod returns the modulus r of the division x / y for y > 0,
+// with x = y*q + r and 0 <= r < y.
+// If y == 0, a division-by-zero run-time error occurs.
+//
+func (x Natural) Mod(y Natural) Natural {
+	_, r := divmod(unpack(x), unpack(y))
+	return pack(r)
+}
+
+
+// DivMod returns the pair (x.Div(y), x.Mod(y)) for y > 0.
+// If y == 0, a division-by-zero run-time error occurs.
+//
+func (x Natural) DivMod(y Natural) (Natural, Natural) {
+	q, r := divmod(unpack(x), unpack(y))
+	return pack(q), pack(r)
+}
+
+
+func shl(z, x Natural, s uint) digit {
+	assert(s <= _W)
+	n := len(x)
+	c := digit(0)
+	for i := 0; i < n; i++ {
+		c, z[i] = x[i]>>(_W-s), x[i]<<s&_M|c
+	}
+	return c
+}
+
+
+// Shl implements ``shift left'' x << s. It returns x * 2^s.
+//
+func (x Natural) Shl(s uint) Natural {
+	n := uint(len(x))
+	m := n + s/_W
+	z := make(Natural, m+1)
+
+	z[m] = shl(z[m-n:m], x, s%_W)
+
+	return normalize(z)
+}
+
+
+func shr(z, x Natural, s uint) digit {
+	assert(s <= _W)
+	n := len(x)
+	c := digit(0)
+	for i := n - 1; i >= 0; i-- {
+		c, z[i] = x[i]<<(_W-s)&_M, x[i]>>s|c
+	}
+	return c
+}
+
+
+// Shr implements ``shift right'' x >> s. It returns x / 2^s.
+//
+func (x Natural) Shr(s uint) Natural {
+	n := uint(len(x))
+	m := n - s/_W
+	if m > n { // check for underflow
+		m = 0
+	}
+	z := make(Natural, m)
+
+	shr(z, x[n-m:n], s%_W)
+
+	return normalize(z)
+}
+
+
+// And returns the ``bitwise and'' x & y for the 2's-complement representation of x and y.
+//
+func (x Natural) And(y Natural) Natural {
+	n := len(x)
+	m := len(y)
+	if n < m {
+		return y.And(x)
+	}
+
+	z := make(Natural, m)
+	for i := 0; i < m; i++ {
+		z[i] = x[i] & y[i]
+	}
+	// upper bits are 0
+
+	return normalize(z)
+}
+
+
+func copy(z, x Natural) {
+	for i, e := range x {
+		z[i] = e
+	}
+}
+
+
+// AndNot returns the ``bitwise clear'' x &^ y for the 2's-complement representation of x and y.
+//
+func (x Natural) AndNot(y Natural) Natural {
+	n := len(x)
+	m := len(y)
+	if n < m {
+		m = n
+	}
+
+	z := make(Natural, n)
+	for i := 0; i < m; i++ {
+		z[i] = x[i] &^ y[i]
+	}
+	copy(z[m:n], x[m:n])
+
+	return normalize(z)
+}
+
+
+// Or returns the ``bitwise or'' x | y for the 2's-complement representation of x and y.
+//
+func (x Natural) Or(y Natural) Natural {
+	n := len(x)
+	m := len(y)
+	if n < m {
+		return y.Or(x)
+	}
+
+	z := make(Natural, n)
+	for i := 0; i < m; i++ {
+		z[i] = x[i] | y[i]
+	}
+	copy(z[m:n], x[m:n])
+
+	return z
+}
+
+
+// Xor returns the ``bitwise exclusive or'' x ^ y for the 2's-complement representation of x and y.
+//
+func (x Natural) Xor(y Natural) Natural {
+	n := len(x)
+	m := len(y)
+	if n < m {
+		return y.Xor(x)
+	}
+
+	z := make(Natural, n)
+	for i := 0; i < m; i++ {
+		z[i] = x[i] ^ y[i]
+	}
+	copy(z[m:n], x[m:n])
+
+	return normalize(z)
+}
+
+
+// Cmp compares x and y. The result is an int value
+//
+//   <  0 if x <  y
+//   == 0 if x == y
+//   >  0 if x >  y
+//
+func (x Natural) Cmp(y Natural) int {
+	n := len(x)
+	m := len(y)
+
+	if n != m || n == 0 {
+		return n - m
+	}
+
+	i := n - 1
+	for i > 0 && x[i] == y[i] {
+		i--
+	}
+
+	d := 0
+	switch {
+	case x[i] < y[i]:
+		d = -1
+	case x[i] > y[i]:
+		d = 1
+	}
+
+	return d
+}
+
+
+// log2 computes the binary logarithm of x for x > 0.
+// The result is the integer n for which 2^n <= x < 2^(n+1).
+// If x == 0 a run-time error occurs.
+//
+func log2(x uint64) uint {
+	assert(x > 0)
+	n := uint(0)
+	for x > 0 {
+		x >>= 1
+		n++
+	}
+	return n - 1
+}
+
+
+// Log2 computes the binary logarithm of x for x > 0.
+// The result is the integer n for which 2^n <= x < 2^(n+1).
+// If x == 0 a run-time error occurs.
+//
+func (x Natural) Log2() uint {
+	n := len(x)
+	if n > 0 {
+		return (uint(n)-1)*_W + log2(uint64(x[n-1]))
+	}
+	panic("Log2(0)")
+}
+
+
+// Computes x = x div d in place (modifies x) for small d's.
+// Returns updated x and x mod d.
+//
+func divmod1(x Natural, d digit) (Natural, digit) {
+	assert(0 < d && isSmall(d-1))
+
+	c := digit(0)
+	for i := len(x) - 1; i >= 0; i-- {
+		t := c<<_W + x[i]
+		c, x[i] = t%d, t/d
+	}
+
+	return normalize(x), c
+}
+
+
+// ToString converts x to a string for a given base, with 2 <= base <= 16.
+//
+func (x Natural) ToString(base uint) string {
+	if len(x) == 0 {
+		return "0"
+	}
+
+	// allocate buffer for conversion
+	assert(2 <= base && base <= 16)
+	n := (x.Log2()+1)/log2(uint64(base)) + 1 // +1: round up
+	s := make([]byte, n)
+
+	// don't destroy x
+	t := make(Natural, len(x))
+	copy(t, x)
+
+	// convert
+	i := n
+	for !t.IsZero() {
+		i--
+		var d digit
+		t, d = divmod1(t, digit(base))
+		s[i] = "0123456789abcdef"[d]
+	}
+
+	return string(s[i:n])
+}
+
+
+// String converts x to its decimal string representation.
+// x.String() is the same as x.ToString(10).
+//
+func (x Natural) String() string { return x.ToString(10) }
+
+
+func fmtbase(c int) uint {
+	switch c {
+	case 'b':
+		return 2
+	case 'o':
+		return 8
+	case 'x':
+		return 16
+	}
+	return 10
+}
+
+
+// Format is a support routine for fmt.Formatter. It accepts
+// the formats 'b' (binary), 'o' (octal), and 'x' (hexadecimal).
+//
+func (x Natural) Format(h fmt.State, c int) { fmt.Fprintf(h, "%s", x.ToString(fmtbase(c))) }
+
+
+func hexvalue(ch byte) uint {
+	d := uint(1 << logH)
+	switch {
+	case '0' <= ch && ch <= '9':
+		d = uint(ch - '0')
+	case 'a' <= ch && ch <= 'f':
+		d = uint(ch-'a') + 10
+	case 'A' <= ch && ch <= 'F':
+		d = uint(ch-'A') + 10
+	}
+	return d
+}
+
+
+// NatFromString returns the natural number corresponding to the
+// longest possible prefix of s representing a natural number in a
+// given conversion base, the actual conversion base used, and the
+// prefix length. The syntax of natural numbers follows the syntax
+// of unsigned integer literals in Go.
+//
+// If the base argument is 0, the string prefix determines the actual
+// conversion base. A prefix of ``0x'' or ``0X'' selects base 16; the
+// ``0'' prefix selects base 8. Otherwise the selected base is 10.
+//
+func NatFromString(s string, base uint) (Natural, uint, int) {
+	// determine base if necessary
+	i, n := 0, len(s)
+	if base == 0 {
+		base = 10
+		if n > 0 && s[0] == '0' {
+			if n > 1 && (s[1] == 'x' || s[1] == 'X') {
+				base, i = 16, 2
+			} else {
+				base, i = 8, 1
+			}
+		}
+	}
+
+	// convert string
+	assert(2 <= base && base <= 16)
+	x := Nat(0)
+	for ; i < n; i++ {
+		d := hexvalue(s[i])
+		if d < base {
+			x = muladd1(x, digit(base), digit(d))
+		} else {
+			break
+		}
+	}
+
+	return x, base, i
+}
+
+
+// Natural number functions
+
+func pop1(x digit) uint {
+	n := uint(0)
+	for x != 0 {
+		x &= x - 1
+		n++
+	}
+	return n
+}
+
+
+// Pop computes the ``population count'' of (the number of 1 bits in) x.
+//
+func (x Natural) Pop() uint {
+	n := uint(0)
+	for i := len(x) - 1; i >= 0; i-- {
+		n += pop1(x[i])
+	}
+	return n
+}
+
+
+// Pow computes x to the power of n.
+//
+func (xp Natural) Pow(n uint) Natural {
+	z := Nat(1)
+	x := xp
+	for n > 0 {
+		// z * x^n == x^n0
+		if n&1 == 1 {
+			z = z.Mul(x)
+		}
+		x, n = x.Mul(x), n/2
+	}
+	return z
+}
+
+
+// MulRange computes the product of all the unsigned integers
+// in the range [a, b] inclusively.
+//
+func MulRange(a, b uint) Natural {
+	switch {
+	case a > b:
+		return Nat(1)
+	case a == b:
+		return Nat(uint64(a))
+	case a+1 == b:
+		return Nat(uint64(a)).Mul(Nat(uint64(b)))
+	}
+	m := (a + b) >> 1
+	assert(a <= m && m < b)
+	return MulRange(a, m).Mul(MulRange(m+1, b))
+}
+
+
+// Fact computes the factorial of n (Fact(n) == MulRange(2, n)).
+//
+func Fact(n uint) Natural {
+	// Using MulRange() instead of the basic for-loop
+	// lead to faster factorial computation.
+	return MulRange(2, n)
+}
+
+
+// Binomial computes the binomial coefficient of (n, k).
+//
+func Binomial(n, k uint) Natural { return MulRange(n-k+1, n).Div(MulRange(1, k)) }
+
+
+// Gcd computes the gcd of x and y.
+//
+func (x Natural) Gcd(y Natural) Natural {
+	// Euclidean algorithm.
+	a, b := x, y
+	for !b.IsZero() {
+		a, b = b, a.Mod(b)
+	}
+	return a
+}
diff --git a/src/pkg/exp/bignum/bignum_test.go b/src/pkg/exp/bignum/bignum_test.go
new file mode 100644
index 000000000..8db93aa96
--- /dev/null
+++ b/src/pkg/exp/bignum/bignum_test.go
@@ -0,0 +1,681 @@
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package bignum
+
+import (
+	"fmt"
+	"testing"
+)
+
+const (
+	sa = "991"
+	sb = "2432902008176640000" // 20!
+	sc = "933262154439441526816992388562667004907159682643816214685929" +
+		"638952175999932299156089414639761565182862536979208272237582" +
+		"51185210916864000000000000000000000000" // 100!
+	sp = "170141183460469231731687303715884105727" // prime
+)
+
+func natFromString(s string, base uint, slen *int) Natural {
+	x, _, len := NatFromString(s, base)
+	if slen != nil {
+		*slen = len
+	}
+	return x
+}
+
+
+func intFromString(s string, base uint, slen *int) *Integer {
+	x, _, len := IntFromString(s, base)
+	if slen != nil {
+		*slen = len
+	}
+	return x
+}
+
+
+func ratFromString(s string, base uint, slen *int) *Rational {
+	x, _, len := RatFromString(s, base)
+	if slen != nil {
+		*slen = len
+	}
+	return x
+}
+
+
+var (
+	nat_zero = Nat(0)
+	nat_one  = Nat(1)
+	nat_two  = Nat(2)
+	a        = natFromString(sa, 10, nil)
+	b        = natFromString(sb, 10, nil)
+	c        = natFromString(sc, 10, nil)
+	p        = natFromString(sp, 10, nil)
+	int_zero = Int(0)
+	int_one  = Int(1)
+	int_two  = Int(2)
+	ip       = intFromString(sp, 10, nil)
+	rat_zero = Rat(0, 1)
+	rat_half = Rat(1, 2)
+	rat_one  = Rat(1, 1)
+	rat_two  = Rat(2, 1)
+)
+
+
+var test_msg string
+var tester *testing.T
+
+func test(n uint, b bool) {
+	if !b {
+		tester.Fatalf("TEST failed: %s (%d)", test_msg, n)
+	}
+}
+
+
+func nat_eq(n uint, x, y Natural) {
+	if x.Cmp(y) != 0 {
+		tester.Fatalf("TEST failed: %s (%d)\nx = %v\ny = %v", test_msg, n, &x, &y)
+	}
+}
+
+
+func int_eq(n uint, x, y *Integer) {
+	if x.Cmp(y) != 0 {
+		tester.Fatalf("TEST failed: %s (%d)\nx = %v\ny = %v", test_msg, n, x, y)
+	}
+}
+
+
+func rat_eq(n uint, x, y *Rational) {
+	if x.Cmp(y) != 0 {
+		tester.Fatalf("TEST failed: %s (%d)\nx = %v\ny = %v", test_msg, n, x, y)
+	}
+}
+
+
+func TestNatConv(t *testing.T) {
+	tester = t
+	test_msg = "NatConvA"
+	type entry1 struct {
+		x uint64
+		s string
+	}
+	tab := []entry1{
+		entry1{0, "0"},
+		entry1{255, "255"},
+		entry1{65535, "65535"},
+		entry1{4294967295, "4294967295"},
+		entry1{18446744073709551615, "18446744073709551615"},
+	}
+	for i, e := range tab {
+		test(100+uint(i), Nat(e.x).String() == e.s)
+		test(200+uint(i), natFromString(e.s, 0, nil).Value() == e.x)
+	}
+
+	test_msg = "NatConvB"
+	for i := uint(0); i < 100; i++ {
+		test(i, Nat(uint64(i)).String() == fmt.Sprintf("%d", i))
+	}
+
+	test_msg = "NatConvC"
+	z := uint64(7)
+	for i := uint(0); i <= 64; i++ {
+		test(i, Nat(z).Value() == z)
+		z <<= 1
+	}
+
+	test_msg = "NatConvD"
+	nat_eq(0, a, Nat(991))
+	nat_eq(1, b, Fact(20))
+	nat_eq(2, c, Fact(100))
+	test(3, a.String() == sa)
+	test(4, b.String() == sb)
+	test(5, c.String() == sc)
+
+	test_msg = "NatConvE"
+	var slen int
+	nat_eq(10, natFromString("0", 0, nil), nat_zero)
+	nat_eq(11, natFromString("123", 0, nil), Nat(123))
+	nat_eq(12, natFromString("077", 0, nil), Nat(7*8+7))
+	nat_eq(13, natFromString("0x1f", 0, nil), Nat(1*16+15))
+	nat_eq(14, natFromString("0x1fg", 0, &slen), Nat(1*16+15))
+	test(4, slen == 4)
+
+	test_msg = "NatConvF"
+	tmp := c.Mul(c)
+	for base := uint(2); base <= 16; base++ {
+		nat_eq(base, natFromString(tmp.ToString(base), base, nil), tmp)
+	}
+
+	test_msg = "NatConvG"
+	x := Nat(100)
+	y, _, _ := NatFromString(fmt.Sprintf("%b", &x), 2)
+	nat_eq(100, y, x)
+}
+
+
+func abs(x int64) uint64 {
+	if x < 0 {
+		x = -x
+	}
+	return uint64(x)
+}
+
+
+func TestIntConv(t *testing.T) {
+	tester = t
+	test_msg = "IntConvA"
+	type entry2 struct {
+		x int64
+		s string
+	}
+	tab := []entry2{
+		entry2{0, "0"},
+		entry2{-128, "-128"},
+		entry2{127, "127"},
+		entry2{-32768, "-32768"},
+		entry2{32767, "32767"},
+		entry2{-2147483648, "-2147483648"},
+		entry2{2147483647, "2147483647"},
+		entry2{-9223372036854775808, "-9223372036854775808"},
+		entry2{9223372036854775807, "9223372036854775807"},
+	}
+	for i, e := range tab {
+		test(100+uint(i), Int(e.x).String() == e.s)
+		test(200+uint(i), intFromString(e.s, 0, nil).Value() == e.x)
+		test(300+uint(i), Int(e.x).Abs().Value() == abs(e.x))
+	}
+
+	test_msg = "IntConvB"
+	var slen int
+	int_eq(0, intFromString("0", 0, nil), int_zero)
+	int_eq(1, intFromString("-0", 0, nil), int_zero)
+	int_eq(2, intFromString("123", 0, nil), Int(123))
+	int_eq(3, intFromString("-123", 0, nil), Int(-123))
+	int_eq(4, intFromString("077", 0, nil), Int(7*8+7))
+	int_eq(5, intFromString("-077", 0, nil), Int(-(7*8 + 7)))
+	int_eq(6, intFromString("0x1f", 0, nil), Int(1*16+15))
+	int_eq(7, intFromString("-0x1f", 0, &slen), Int(-(1*16 + 15)))
+	test(7, slen == 5)
+	int_eq(8, intFromString("+0x1f", 0, &slen), Int(+(1*16 + 15)))
+	test(8, slen == 5)
+	int_eq(9, intFromString("0x1fg", 0, &slen), Int(1*16+15))
+	test(9, slen == 4)
+	int_eq(10, intFromString("-0x1fg", 0, &slen), Int(-(1*16 + 15)))
+	test(10, slen == 5)
+}
+
+
+func TestRatConv(t *testing.T) {
+	tester = t
+	test_msg = "RatConv"
+	var slen int
+	rat_eq(0, ratFromString("0", 0, nil), rat_zero)
+	rat_eq(1, ratFromString("0/1", 0, nil), rat_zero)
+	rat_eq(2, ratFromString("0/01", 0, nil), rat_zero)
+	rat_eq(3, ratFromString("0x14/10", 0, &slen), rat_two)
+	test(4, slen == 7)
+	rat_eq(5, ratFromString("0.", 0, nil), rat_zero)
+	rat_eq(6, ratFromString("0.001f", 10, nil), Rat(1, 1000))
+	rat_eq(7, ratFromString(".1", 0, nil), Rat(1, 10))
+	rat_eq(8, ratFromString("10101.0101", 2, nil), Rat(0x155, 1<<4))
+	rat_eq(9, ratFromString("-0003.145926", 10, &slen), Rat(-3145926, 1000000))
+	test(10, slen == 12)
+	rat_eq(11, ratFromString("1e2", 0, nil), Rat(100, 1))
+	rat_eq(12, ratFromString("1e-2", 0, nil), Rat(1, 100))
+	rat_eq(13, ratFromString("1.1e2", 0, nil), Rat(110, 1))
+	rat_eq(14, ratFromString(".1e2x", 0, &slen), Rat(10, 1))
+	test(15, slen == 4)
+}
+
+
+func add(x, y Natural) Natural {
+	z1 := x.Add(y)
+	z2 := y.Add(x)
+	if z1.Cmp(z2) != 0 {
+		tester.Fatalf("addition not symmetric:\n\tx = %v\n\ty = %t", x, y)
+	}
+	return z1
+}
+
+
+func sum(n uint64, scale Natural) Natural {
+	s := nat_zero
+	for ; n > 0; n-- {
+		s = add(s, Nat(n).Mul(scale))
+	}
+	return s
+}
+
+
+func TestNatAdd(t *testing.T) {
+	tester = t
+	test_msg = "NatAddA"
+	nat_eq(0, add(nat_zero, nat_zero), nat_zero)
+	nat_eq(1, add(nat_zero, c), c)
+
+	test_msg = "NatAddB"
+	for i := uint64(0); i < 100; i++ {
+		t := Nat(i)
+		nat_eq(uint(i), sum(i, c), t.Mul(t).Add(t).Shr(1).Mul(c))
+	}
+}
+
+
+func mul(x, y Natural) Natural {
+	z1 := x.Mul(y)
+	z2 := y.Mul(x)
+	if z1.Cmp(z2) != 0 {
+		tester.Fatalf("multiplication not symmetric:\n\tx = %v\n\ty = %t", x, y)
+	}
+	if !x.IsZero() && z1.Div(x).Cmp(y) != 0 {
+		tester.Fatalf("multiplication/division not inverse (A):\n\tx = %v\n\ty = %t", x, y)
+	}
+	if !y.IsZero() && z1.Div(y).Cmp(x) != 0 {
+		tester.Fatalf("multiplication/division not inverse (B):\n\tx = %v\n\ty = %t", x, y)
+	}
+	return z1
+}
+
+
+func TestNatSub(t *testing.T) {
+	tester = t
+	test_msg = "NatSubA"
+	nat_eq(0, nat_zero.Sub(nat_zero), nat_zero)
+	nat_eq(1, c.Sub(nat_zero), c)
+
+	test_msg = "NatSubB"
+	for i := uint64(0); i < 100; i++ {
+		t := sum(i, c)
+		for j := uint64(0); j <= i; j++ {
+			t = t.Sub(mul(Nat(j), c))
+		}
+		nat_eq(uint(i), t, nat_zero)
+	}
+}
+
+
+func TestNatMul(t *testing.T) {
+	tester = t
+	test_msg = "NatMulA"
+	nat_eq(0, mul(c, nat_zero), nat_zero)
+	nat_eq(1, mul(c, nat_one), c)
+
+	test_msg = "NatMulB"
+	nat_eq(0, b.Mul(MulRange(0, 100)), nat_zero)
+	nat_eq(1, b.Mul(MulRange(21, 100)), c)
+
+	test_msg = "NatMulC"
+	const n = 100
+	p := b.Mul(c).Shl(n)
+	for i := uint(0); i < n; i++ {
+		nat_eq(i, mul(b.Shl(i), c.Shl(n-i)), p)
+	}
+}
+
+
+func TestNatDiv(t *testing.T) {
+	tester = t
+	test_msg = "NatDivA"
+	nat_eq(0, c.Div(nat_one), c)
+	nat_eq(1, c.Div(Nat(100)), Fact(99))
+	nat_eq(2, b.Div(c), nat_zero)
+	nat_eq(4, nat_one.Shl(100).Div(nat_one.Shl(90)), nat_one.Shl(10))
+	nat_eq(5, c.Div(b), MulRange(21, 100))
+
+	test_msg = "NatDivB"
+	const n = 100
+	p := Fact(n)
+	for i := uint(0); i < n; i++ {
+		nat_eq(100+i, p.Div(MulRange(1, i)), MulRange(i+1, n))
+	}
+
+	// a specific test case that exposed a bug in package big
+	test_msg = "NatDivC"
+	x := natFromString("69720375229712477164533808935312303556800", 10, nil)
+	y := natFromString("3099044504245996706400", 10, nil)
+	q := natFromString("22497377864108980962", 10, nil)
+	r := natFromString("0", 10, nil)
+	qc, rc := x.DivMod(y)
+	nat_eq(0, q, qc)
+	nat_eq(1, r, rc)
+}
+
+
+func TestIntQuoRem(t *testing.T) {
+	tester = t
+	test_msg = "IntQuoRem"
+	type T struct {
+		x, y, q, r int64
+	}
+	a := []T{
+		T{+8, +3, +2, +2},
+		T{+8, -3, -2, +2},
+		T{-8, +3, -2, -2},
+		T{-8, -3, +2, -2},
+		T{+1, +2, 0, +1},
+		T{+1, -2, 0, +1},
+		T{-1, +2, 0, -1},
+		T{-1, -2, 0, -1},
+	}
+	for i := uint(0); i < uint(len(a)); i++ {
+		e := &a[i]
+		x, y := Int(e.x).Mul(ip), Int(e.y).Mul(ip)
+		q, r := Int(e.q), Int(e.r).Mul(ip)
+		qq, rr := x.QuoRem(y)
+		int_eq(4*i+0, x.Quo(y), q)
+		int_eq(4*i+1, x.Rem(y), r)
+		int_eq(4*i+2, qq, q)
+		int_eq(4*i+3, rr, r)
+	}
+}
+
+
+func TestIntDivMod(t *testing.T) {
+	tester = t
+	test_msg = "IntDivMod"
+	type T struct {
+		x, y, q, r int64
+	}
+	a := []T{
+		T{+8, +3, +2, +2},
+		T{+8, -3, -2, +2},
+		T{-8, +3, -3, +1},
+		T{-8, -3, +3, +1},
+		T{+1, +2, 0, +1},
+		T{+1, -2, 0, +1},
+		T{-1, +2, -1, +1},
+		T{-1, -2, +1, +1},
+	}
+	for i := uint(0); i < uint(len(a)); i++ {
+		e := &a[i]
+		x, y := Int(e.x).Mul(ip), Int(e.y).Mul(ip)
+		q, r := Int(e.q), Int(e.r).Mul(ip)
+		qq, rr := x.DivMod(y)
+		int_eq(4*i+0, x.Div(y), q)
+		int_eq(4*i+1, x.Mod(y), r)
+		int_eq(4*i+2, qq, q)
+		int_eq(4*i+3, rr, r)
+	}
+}
+
+
+func TestNatMod(t *testing.T) {
+	tester = t
+	test_msg = "NatModA"
+	for i := uint(0); ; i++ {
+		d := nat_one.Shl(i)
+		if d.Cmp(c) < 0 {
+			nat_eq(i, c.Add(d).Mod(c), d)
+		} else {
+			nat_eq(i, c.Add(d).Div(c), nat_two)
+			nat_eq(i, c.Add(d).Mod(c), d.Sub(c))
+			break
+		}
+	}
+}
+
+
+func TestNatShift(t *testing.T) {
+	tester = t
+	test_msg = "NatShift1L"
+	test(0, b.Shl(0).Cmp(b) == 0)
+	test(1, c.Shl(1).Cmp(c) > 0)
+
+	test_msg = "NatShift1R"
+	test(3, b.Shr(0).Cmp(b) == 0)
+	test(4, c.Shr(1).Cmp(c) < 0)
+
+	test_msg = "NatShift2"
+	for i := uint(0); i < 100; i++ {
+		test(i, c.Shl(i).Shr(i).Cmp(c) == 0)
+	}
+
+	test_msg = "NatShift3L"
+	{
+		const m = 3
+		p := b
+		f := Nat(1 << m)
+		for i := uint(0); i < 100; i++ {
+			nat_eq(i, b.Shl(i*m), p)
+			p = mul(p, f)
+		}
+	}
+
+	test_msg = "NatShift3R"
+	{
+		p := c
+		for i := uint(0); !p.IsZero(); i++ {
+			nat_eq(i, c.Shr(i), p)
+			p = p.Shr(1)
+		}
+	}
+}
+
+
+func TestIntShift(t *testing.T) {
+	tester = t
+	test_msg = "IntShift1L"
+	test(0, ip.Shl(0).Cmp(ip) == 0)
+	test(1, ip.Shl(1).Cmp(ip) > 0)
+
+	test_msg = "IntShift1R"
+	test(0, ip.Shr(0).Cmp(ip) == 0)
+	test(1, ip.Shr(1).Cmp(ip) < 0)
+
+	test_msg = "IntShift2"
+	for i := uint(0); i < 100; i++ {
+		test(i, ip.Shl(i).Shr(i).Cmp(ip) == 0)
+	}
+
+	test_msg = "IntShift3L"
+	{
+		const m = 3
+		p := ip
+		f := Int(1 << m)
+		for i := uint(0); i < 100; i++ {
+			int_eq(i, ip.Shl(i*m), p)
+			p = p.Mul(f)
+		}
+	}
+
+	test_msg = "IntShift3R"
+	{
+		p := ip
+		for i := uint(0); p.IsPos(); i++ {
+			int_eq(i, ip.Shr(i), p)
+			p = p.Shr(1)
+		}
+	}
+
+	test_msg = "IntShift4R"
+	int_eq(0, Int(-43).Shr(1), Int(-43>>1))
+	int_eq(0, Int(-1024).Shr(100), Int(-1))
+	int_eq(1, ip.Neg().Shr(10), ip.Neg().Div(Int(1).Shl(10)))
+}
+
+
+func TestNatBitOps(t *testing.T) {
+	tester = t
+
+	x := uint64(0xf08e6f56bd8c3941)
+	y := uint64(0x3984ef67834bc)
+
+	bx := Nat(x)
+	by := Nat(y)
+
+	test_msg = "NatAnd"
+	bz := Nat(x & y)
+	for i := uint(0); i < 100; i++ {
+		nat_eq(i, bx.Shl(i).And(by.Shl(i)), bz.Shl(i))
+	}
+
+	test_msg = "NatAndNot"
+	bz = Nat(x &^ y)
+	for i := uint(0); i < 100; i++ {
+		nat_eq(i, bx.Shl(i).AndNot(by.Shl(i)), bz.Shl(i))
+	}
+
+	test_msg = "NatOr"
+	bz = Nat(x | y)
+	for i := uint(0); i < 100; i++ {
+		nat_eq(i, bx.Shl(i).Or(by.Shl(i)), bz.Shl(i))
+	}
+
+	test_msg = "NatXor"
+	bz = Nat(x ^ y)
+	for i := uint(0); i < 100; i++ {
+		nat_eq(i, bx.Shl(i).Xor(by.Shl(i)), bz.Shl(i))
+	}
+}
+
+
+func TestIntBitOps1(t *testing.T) {
+	tester = t
+	test_msg = "IntBitOps1"
+	type T struct {
+		x, y int64
+	}
+	a := []T{
+		T{+7, +3},
+		T{+7, -3},
+		T{-7, +3},
+		T{-7, -3},
+	}
+	for i := uint(0); i < uint(len(a)); i++ {
+		e := &a[i]
+		int_eq(4*i+0, Int(e.x).And(Int(e.y)), Int(e.x&e.y))
+		int_eq(4*i+1, Int(e.x).AndNot(Int(e.y)), Int(e.x&^e.y))
+		int_eq(4*i+2, Int(e.x).Or(Int(e.y)), Int(e.x|e.y))
+		int_eq(4*i+3, Int(e.x).Xor(Int(e.y)), Int(e.x^e.y))
+	}
+}
+
+
+func TestIntBitOps2(t *testing.T) {
+	tester = t
+
+	test_msg = "IntNot"
+	int_eq(0, Int(-2).Not(), Int(1))
+	int_eq(0, Int(-1).Not(), Int(0))
+	int_eq(0, Int(0).Not(), Int(-1))
+	int_eq(0, Int(1).Not(), Int(-2))
+	int_eq(0, Int(2).Not(), Int(-3))
+
+	test_msg = "IntAnd"
+	for x := int64(-15); x < 5; x++ {
+		bx := Int(x)
+		for y := int64(-5); y < 15; y++ {
+			by := Int(y)
+			for i := uint(50); i < 70; i++ { // shift across 64bit boundary
+				int_eq(i, bx.Shl(i).And(by.Shl(i)), Int(x&y).Shl(i))
+			}
+		}
+	}
+
+	test_msg = "IntAndNot"
+	for x := int64(-15); x < 5; x++ {
+		bx := Int(x)
+		for y := int64(-5); y < 15; y++ {
+			by := Int(y)
+			for i := uint(50); i < 70; i++ { // shift across 64bit boundary
+				int_eq(2*i+0, bx.Shl(i).AndNot(by.Shl(i)), Int(x&^y).Shl(i))
+				int_eq(2*i+1, bx.Shl(i).And(by.Shl(i).Not()), Int(x&^y).Shl(i))
+			}
+		}
+	}
+
+	test_msg = "IntOr"
+	for x := int64(-15); x < 5; x++ {
+		bx := Int(x)
+		for y := int64(-5); y < 15; y++ {
+			by := Int(y)
+			for i := uint(50); i < 70; i++ { // shift across 64bit boundary
+				int_eq(i, bx.Shl(i).Or(by.Shl(i)), Int(x|y).Shl(i))
+			}
+		}
+	}
+
+	test_msg = "IntXor"
+	for x := int64(-15); x < 5; x++ {
+		bx := Int(x)
+		for y := int64(-5); y < 15; y++ {
+			by := Int(y)
+			for i := uint(50); i < 70; i++ { // shift across 64bit boundary
+				int_eq(i, bx.Shl(i).Xor(by.Shl(i)), Int(x^y).Shl(i))
+			}
+		}
+	}
+}
+
+
+func TestNatCmp(t *testing.T) {
+	tester = t
+	test_msg = "NatCmp"
+	test(0, a.Cmp(a) == 0)
+	test(1, a.Cmp(b) < 0)
+	test(2, b.Cmp(a) > 0)
+	test(3, a.Cmp(c) < 0)
+	d := c.Add(b)
+	test(4, c.Cmp(d) < 0)
+	test(5, d.Cmp(c) > 0)
+}
+
+
+func TestNatLog2(t *testing.T) {
+	tester = t
+	test_msg = "NatLog2A"
+	test(0, nat_one.Log2() == 0)
+	test(1, nat_two.Log2() == 1)
+	test(2, Nat(3).Log2() == 1)
+	test(3, Nat(4).Log2() == 2)
+
+	test_msg = "NatLog2B"
+	for i := uint(0); i < 100; i++ {
+		test(i, nat_one.Shl(i).Log2() == i)
+	}
+}
+
+
+func TestNatGcd(t *testing.T) {
+	tester = t
+	test_msg = "NatGcdA"
+	f := Nat(99991)
+	nat_eq(0, b.Mul(f).Gcd(c.Mul(f)), MulRange(1, 20).Mul(f))
+}
+
+
+func TestNatPow(t *testing.T) {
+	tester = t
+	test_msg = "NatPowA"
+	nat_eq(0, nat_two.Pow(0), nat_one)
+
+	test_msg = "NatPowB"
+	for i := uint(0); i < 100; i++ {
+		nat_eq(i, nat_two.Pow(i), nat_one.Shl(i))
+	}
+}
+
+
+func TestNatPop(t *testing.T) {
+	tester = t
+	test_msg = "NatPopA"
+	test(0, nat_zero.Pop() == 0)
+	test(1, nat_one.Pop() == 1)
+	test(2, Nat(10).Pop() == 2)
+	test(3, Nat(30).Pop() == 4)
+	test(4, Nat(0x1248f).Shl(33).Pop() == 8)
+
+	test_msg = "NatPopB"
+	for i := uint(0); i < 100; i++ {
+		test(i, nat_one.Shl(i).Sub(nat_one).Pop() == i)
+	}
+}
+
+
+func TestIssue571(t *testing.T) {
+	const min_float = "4.940656458412465441765687928682213723651e-324"
+	RatFromString(min_float, 10) // this must not crash
+}
diff --git a/src/pkg/exp/bignum/integer.go b/src/pkg/exp/bignum/integer.go
new file mode 100644
index 000000000..a8d26829d
--- /dev/null
+++ b/src/pkg/exp/bignum/integer.go
@@ -0,0 +1,520 @@
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Integer numbers
+//
+// Integers are normalized if the mantissa is normalized and the sign is
+// false for mant == 0. Use MakeInt to create normalized Integers.
+
+package bignum
+
+import (
+	"fmt"
+)
+
+// TODO(gri) Complete the set of in-place operations.
+
+// Integer represents a signed integer value of arbitrary precision.
+//
+type Integer struct {
+	sign bool
+	mant Natural
+}
+
+
+// MakeInt makes an integer given a sign and a mantissa.
+// The number is positive (>= 0) if sign is false or the
+// mantissa is zero; it is negative otherwise.
+//
+func MakeInt(sign bool, mant Natural) *Integer {
+	if mant.IsZero() {
+		sign = false // normalize
+	}
+	return &Integer{sign, mant}
+}
+
+
+// Int creates a small integer with value x.
+//
+func Int(x int64) *Integer {
+	var ux uint64
+	if x < 0 {
+		// For the most negative x, -x == x, and
+		// the bit pattern has the correct value.
+		ux = uint64(-x)
+	} else {
+		ux = uint64(x)
+	}
+	return MakeInt(x < 0, Nat(ux))
+}
+
+
+// Value returns the value of x, if x fits into an int64;
+// otherwise the result is undefined.
+//
+func (x *Integer) Value() int64 {
+	z := int64(x.mant.Value())
+	if x.sign {
+		z = -z
+	}
+	return z
+}
+
+
+// Abs returns the absolute value of x.
+//
+func (x *Integer) Abs() Natural { return x.mant }
+
+
+// Predicates
+
+// IsEven returns true iff x is divisible by 2.
+//
+func (x *Integer) IsEven() bool { return x.mant.IsEven() }
+
+
+// IsOdd returns true iff x is not divisible by 2.
+//
+func (x *Integer) IsOdd() bool { return x.mant.IsOdd() }
+
+
+// IsZero returns true iff x == 0.
+//
+func (x *Integer) IsZero() bool { return x.mant.IsZero() }
+
+
+// IsNeg returns true iff x < 0.
+//
+func (x *Integer) IsNeg() bool { return x.sign && !x.mant.IsZero() }
+
+
+// IsPos returns true iff x >= 0.
+//
+func (x *Integer) IsPos() bool { return !x.sign && !x.mant.IsZero() }
+
+
+// Operations
+
+// Neg returns the negated value of x.
+//
+func (x *Integer) Neg() *Integer { return MakeInt(!x.sign, x.mant) }
+
+
+// Iadd sets z to the sum x + y.
+// z must exist and may be x or y.
+//
+func Iadd(z, x, y *Integer) {
+	if x.sign == y.sign {
+		// x + y == x + y
+		// (-x) + (-y) == -(x + y)
+		z.sign = x.sign
+		Nadd(&z.mant, x.mant, y.mant)
+	} else {
+		// x + (-y) == x - y == -(y - x)
+		// (-x) + y == y - x == -(x - y)
+		if x.mant.Cmp(y.mant) >= 0 {
+			z.sign = x.sign
+			Nsub(&z.mant, x.mant, y.mant)
+		} else {
+			z.sign = !x.sign
+			Nsub(&z.mant, y.mant, x.mant)
+		}
+	}
+}
+
+
+// Add returns the sum x + y.
+//
+func (x *Integer) Add(y *Integer) *Integer {
+	var z Integer
+	Iadd(&z, x, y)
+	return &z
+}
+
+
+func Isub(z, x, y *Integer) {
+	if x.sign != y.sign {
+		// x - (-y) == x + y
+		// (-x) - y == -(x + y)
+		z.sign = x.sign
+		Nadd(&z.mant, x.mant, y.mant)
+	} else {
+		// x - y == x - y == -(y - x)
+		// (-x) - (-y) == y - x == -(x - y)
+		if x.mant.Cmp(y.mant) >= 0 {
+			z.sign = x.sign
+			Nsub(&z.mant, x.mant, y.mant)
+		} else {
+			z.sign = !x.sign
+			Nsub(&z.mant, y.mant, x.mant)
+		}
+	}
+}
+
+
+// Sub returns the difference x - y.
+//
+func (x *Integer) Sub(y *Integer) *Integer {
+	var z Integer
+	Isub(&z, x, y)
+	return &z
+}
+
+
+// Nscale sets *z to the scaled value (*z) * d.
+//
+func Iscale(z *Integer, d int64) {
+	f := uint64(d)
+	if d < 0 {
+		f = uint64(-d)
+	}
+	z.sign = z.sign != (d < 0)
+	Nscale(&z.mant, f)
+}
+
+
+// Mul1 returns the product x * d.
+//
+func (x *Integer) Mul1(d int64) *Integer {
+	f := uint64(d)
+	if d < 0 {
+		f = uint64(-d)
+	}
+	return MakeInt(x.sign != (d < 0), x.mant.Mul1(f))
+}
+
+
+// Mul returns the product x * y.
+//
+func (x *Integer) Mul(y *Integer) *Integer {
+	// x * y == x * y
+	// x * (-y) == -(x * y)
+	// (-x) * y == -(x * y)
+	// (-x) * (-y) == x * y
+	return MakeInt(x.sign != y.sign, x.mant.Mul(y.mant))
+}
+
+
+// MulNat returns the product x * y, where y is a (unsigned) natural number.
+//
+func (x *Integer) MulNat(y Natural) *Integer {
+	// x * y == x * y
+	// (-x) * y == -(x * y)
+	return MakeInt(x.sign, x.mant.Mul(y))
+}
+
+
+// Quo returns the quotient q = x / y for y != 0.
+// If y == 0, a division-by-zero run-time error occurs.
+//
+// Quo and Rem implement T-division and modulus (like C99):
+//
+//   q = x.Quo(y) = trunc(x/y)  (truncation towards zero)
+//   r = x.Rem(y) = x - y*q
+//
+// (Daan Leijen, ``Division and Modulus for Computer Scientists''.)
+//
+func (x *Integer) Quo(y *Integer) *Integer {
+	// x / y == x / y
+	// x / (-y) == -(x / y)
+	// (-x) / y == -(x / y)
+	// (-x) / (-y) == x / y
+	return MakeInt(x.sign != y.sign, x.mant.Div(y.mant))
+}
+
+
+// Rem returns the remainder r of the division x / y for y != 0,
+// with r = x - y*x.Quo(y). Unless r is zero, its sign corresponds
+// to the sign of x.
+// If y == 0, a division-by-zero run-time error occurs.
+//
+func (x *Integer) Rem(y *Integer) *Integer {
+	// x % y == x % y
+	// x % (-y) == x % y
+	// (-x) % y == -(x % y)
+	// (-x) % (-y) == -(x % y)
+	return MakeInt(x.sign, x.mant.Mod(y.mant))
+}
+
+
+// QuoRem returns the pair (x.Quo(y), x.Rem(y)) for y != 0.
+// If y == 0, a division-by-zero run-time error occurs.
+//
+func (x *Integer) QuoRem(y *Integer) (*Integer, *Integer) {
+	q, r := x.mant.DivMod(y.mant)
+	return MakeInt(x.sign != y.sign, q), MakeInt(x.sign, r)
+}
+
+
+// Div returns the quotient q = x / y for y != 0.
+// If y == 0, a division-by-zero run-time error occurs.
+//
+// Div and Mod implement Euclidian division and modulus:
+//
+//   q = x.Div(y)
+//   r = x.Mod(y) with: 0 <= r < |q| and: x = y*q + r
+//
+// (Raymond T. Boute, ``The Euclidian definition of the functions
+// div and mod''. ACM Transactions on Programming Languages and
+// Systems (TOPLAS), 14(2):127-144, New York, NY, USA, 4/1992.
+// ACM press.)
+//
+func (x *Integer) Div(y *Integer) *Integer {
+	q, r := x.QuoRem(y)
+	if r.IsNeg() {
+		if y.IsPos() {
+			q = q.Sub(Int(1))
+		} else {
+			q = q.Add(Int(1))
+		}
+	}
+	return q
+}
+
+
+// Mod returns the modulus r of the division x / y for y != 0,
+// with r = x - y*x.Div(y). r is always positive.
+// If y == 0, a division-by-zero run-time error occurs.
+//
+func (x *Integer) Mod(y *Integer) *Integer {
+	r := x.Rem(y)
+	if r.IsNeg() {
+		if y.IsPos() {
+			r = r.Add(y)
+		} else {
+			r = r.Sub(y)
+		}
+	}
+	return r
+}
+
+
+// DivMod returns the pair (x.Div(y), x.Mod(y)).
+//
+func (x *Integer) DivMod(y *Integer) (*Integer, *Integer) {
+	q, r := x.QuoRem(y)
+	if r.IsNeg() {
+		if y.IsPos() {
+			q = q.Sub(Int(1))
+			r = r.Add(y)
+		} else {
+			q = q.Add(Int(1))
+			r = r.Sub(y)
+		}
+	}
+	return q, r
+}
+
+
+// Shl implements ``shift left'' x << s. It returns x * 2^s.
+//
+func (x *Integer) Shl(s uint) *Integer { return MakeInt(x.sign, x.mant.Shl(s)) }
+
+
+// The bitwise operations on integers are defined on the 2's-complement
+// representation of integers. From
+//
+//   -x == ^x + 1  (1)  2's complement representation
+//
+// follows:
+//
+//   -(x) == ^(x) + 1
+//   -(-x) == ^(-x) + 1
+//   x-1 == ^(-x)
+//   ^(x-1) == -x  (2)
+//
+// Using (1) and (2), operations on negative integers of the form -x are
+// converted to operations on negated positive integers of the form ~(x-1).
+
+
+// Shr implements ``shift right'' x >> s. It returns x / 2^s.
+//
+func (x *Integer) Shr(s uint) *Integer {
+	if x.sign {
+		// (-x) >> s == ^(x-1) >> s == ^((x-1) >> s) == -(((x-1) >> s) + 1)
+		return MakeInt(true, x.mant.Sub(Nat(1)).Shr(s).Add(Nat(1)))
+	}
+
+	return MakeInt(false, x.mant.Shr(s))
+}
+
+
+// Not returns the ``bitwise not'' ^x for the 2's-complement representation of x.
+func (x *Integer) Not() *Integer {
+	if x.sign {
+		// ^(-x) == ^(^(x-1)) == x-1
+		return MakeInt(false, x.mant.Sub(Nat(1)))
+	}
+
+	// ^x == -x-1 == -(x+1)
+	return MakeInt(true, x.mant.Add(Nat(1)))
+}
+
+
+// And returns the ``bitwise and'' x & y for the 2's-complement representation of x and y.
+//
+func (x *Integer) And(y *Integer) *Integer {
+	if x.sign == y.sign {
+		if x.sign {
+			// (-x) & (-y) == ^(x-1) & ^(y-1) == ^((x-1) | (y-1)) == -(((x-1) | (y-1)) + 1)
+			return MakeInt(true, x.mant.Sub(Nat(1)).Or(y.mant.Sub(Nat(1))).Add(Nat(1)))
+		}
+
+		// x & y == x & y
+		return MakeInt(false, x.mant.And(y.mant))
+	}
+
+	// x.sign != y.sign
+	if x.sign {
+		x, y = y, x // & is symmetric
+	}
+
+	// x & (-y) == x & ^(y-1) == x &^ (y-1)
+	return MakeInt(false, x.mant.AndNot(y.mant.Sub(Nat(1))))
+}
+
+
+// AndNot returns the ``bitwise clear'' x &^ y for the 2's-complement representation of x and y.
+//
+func (x *Integer) AndNot(y *Integer) *Integer {
+	if x.sign == y.sign {
+		if x.sign {
+			// (-x) &^ (-y) == ^(x-1) &^ ^(y-1) == ^(x-1) & (y-1) == (y-1) &^ (x-1)
+			return MakeInt(false, y.mant.Sub(Nat(1)).AndNot(x.mant.Sub(Nat(1))))
+		}
+
+		// x &^ y == x &^ y
+		return MakeInt(false, x.mant.AndNot(y.mant))
+	}
+
+	if x.sign {
+		// (-x) &^ y == ^(x-1) &^ y == ^(x-1) & ^y == ^((x-1) | y) == -(((x-1) | y) + 1)
+		return MakeInt(true, x.mant.Sub(Nat(1)).Or(y.mant).Add(Nat(1)))
+	}
+
+	// x &^ (-y) == x &^ ^(y-1) == x & (y-1)
+	return MakeInt(false, x.mant.And(y.mant.Sub(Nat(1))))
+}
+
+
+// Or returns the ``bitwise or'' x | y for the 2's-complement representation of x and y.
+//
+func (x *Integer) Or(y *Integer) *Integer {
+	if x.sign == y.sign {
+		if x.sign {
+			// (-x) | (-y) == ^(x-1) | ^(y-1) == ^((x-1) & (y-1)) == -(((x-1) & (y-1)) + 1)
+			return MakeInt(true, x.mant.Sub(Nat(1)).And(y.mant.Sub(Nat(1))).Add(Nat(1)))
+		}
+
+		// x | y == x | y
+		return MakeInt(false, x.mant.Or(y.mant))
+	}
+
+	// x.sign != y.sign
+	if x.sign {
+		x, y = y, x // | or symmetric
+	}
+
+	// x | (-y) == x | ^(y-1) == ^((y-1) &^ x) == -(^((y-1) &^ x) + 1)
+	return MakeInt(true, y.mant.Sub(Nat(1)).AndNot(x.mant).Add(Nat(1)))
+}
+
+
+// Xor returns the ``bitwise xor'' x | y for the 2's-complement representation of x and y.
+//
+func (x *Integer) Xor(y *Integer) *Integer {
+	if x.sign == y.sign {
+		if x.sign {
+			// (-x) ^ (-y) == ^(x-1) ^ ^(y-1) == (x-1) ^ (y-1)
+			return MakeInt(false, x.mant.Sub(Nat(1)).Xor(y.mant.Sub(Nat(1))))
+		}
+
+		// x ^ y == x ^ y
+		return MakeInt(false, x.mant.Xor(y.mant))
+	}
+
+	// x.sign != y.sign
+	if x.sign {
+		x, y = y, x // ^ is symmetric
+	}
+
+	// x ^ (-y) == x ^ ^(y-1) == ^(x ^ (y-1)) == -((x ^ (y-1)) + 1)
+	return MakeInt(true, x.mant.Xor(y.mant.Sub(Nat(1))).Add(Nat(1)))
+}
+
+
+// Cmp compares x and y. The result is an int value that is
+//
+//   <  0 if x <  y
+//   == 0 if x == y
+//   >  0 if x >  y
+//
+func (x *Integer) Cmp(y *Integer) int {
+	// x cmp y == x cmp y
+	// x cmp (-y) == x
+	// (-x) cmp y == y
+	// (-x) cmp (-y) == -(x cmp y)
+	var r int
+	switch {
+	case x.sign == y.sign:
+		r = x.mant.Cmp(y.mant)
+		if x.sign {
+			r = -r
+		}
+	case x.sign:
+		r = -1
+	case y.sign:
+		r = 1
+	}
+	return r
+}
+
+
+// ToString converts x to a string for a given base, with 2 <= base <= 16.
+//
+func (x *Integer) ToString(base uint) string {
+	if x.mant.IsZero() {
+		return "0"
+	}
+	var s string
+	if x.sign {
+		s = "-"
+	}
+	return s + x.mant.ToString(base)
+}
+
+
+// String converts x to its decimal string representation.
+// x.String() is the same as x.ToString(10).
+//
+func (x *Integer) String() string { return x.ToString(10) }
+
+
+// Format is a support routine for fmt.Formatter. It accepts
+// the formats 'b' (binary), 'o' (octal), and 'x' (hexadecimal).
+//
+func (x *Integer) Format(h fmt.State, c int) { fmt.Fprintf(h, "%s", x.ToString(fmtbase(c))) }
+
+
+// IntFromString returns the integer corresponding to the
+// longest possible prefix of s representing an integer in a
+// given conversion base, the actual conversion base used, and
+// the prefix length. The syntax of integers follows the syntax
+// of signed integer literals in Go.
+//
+// If the base argument is 0, the string prefix determines the actual
+// conversion base. A prefix of ``0x'' or ``0X'' selects base 16; the
+// ``0'' prefix selects base 8. Otherwise the selected base is 10.
+//
+func IntFromString(s string, base uint) (*Integer, uint, int) {
+	// skip sign, if any
+	i0 := 0
+	if len(s) > 0 && (s[0] == '-' || s[0] == '+') {
+		i0 = 1
+	}
+
+	mant, base, slen := NatFromString(s[i0:], base)
+
+	return MakeInt(i0 > 0 && s[0] == '-', mant), base, i0 + slen
+}
diff --git a/src/pkg/exp/bignum/nrdiv_test.go b/src/pkg/exp/bignum/nrdiv_test.go
new file mode 100644
index 000000000..725b1acea
--- /dev/null
+++ b/src/pkg/exp/bignum/nrdiv_test.go
@@ -0,0 +1,188 @@
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// This file implements Newton-Raphson division and uses
+// it as an additional test case for bignum.
+//
+// Division of x/y is achieved by computing r = 1/y to
+// obtain the quotient q = x*r = x*(1/y) = x/y. The
+// reciprocal r is the solution for f(x) = 1/x - y and
+// the solution is approximated through iteration. The
+// iteration does not require division.
+
+package bignum
+
+import "testing"
+
+
+// An fpNat is a Natural scaled by a power of two
+// (an unsigned floating point representation). The
+// value of an fpNat x is x.m * 2^x.e .
+//
+type fpNat struct {
+	m Natural
+	e int
+}
+
+
+// sub computes x - y.
+func (x fpNat) sub(y fpNat) fpNat {
+	switch d := x.e - y.e; {
+	case d < 0:
+		return fpNat{x.m.Sub(y.m.Shl(uint(-d))), x.e}
+	case d > 0:
+		return fpNat{x.m.Shl(uint(d)).Sub(y.m), y.e}
+	}
+	return fpNat{x.m.Sub(y.m), x.e}
+}
+
+
+// mul2 computes x*2.
+func (x fpNat) mul2() fpNat { return fpNat{x.m, x.e + 1} }
+
+
+// mul computes x*y.
+func (x fpNat) mul(y fpNat) fpNat { return fpNat{x.m.Mul(y.m), x.e + y.e} }
+
+
+// mant computes the (possibly truncated) Natural representation
+// of an fpNat x.
+//
+func (x fpNat) mant() Natural {
+	switch {
+	case x.e > 0:
+		return x.m.Shl(uint(x.e))
+	case x.e < 0:
+		return x.m.Shr(uint(-x.e))
+	}
+	return x.m
+}
+
+
+// nrDivEst computes an estimate of the quotient q = x0/y0 and returns q.
+// q may be too small (usually by 1).
+//
+func nrDivEst(x0, y0 Natural) Natural {
+	if y0.IsZero() {
+		panic("division by zero")
+		return nil
+	}
+	// y0 > 0
+
+	if y0.Cmp(Nat(1)) == 0 {
+		return x0
+	}
+	// y0 > 1
+
+	switch d := x0.Cmp(y0); {
+	case d < 0:
+		return Nat(0)
+	case d == 0:
+		return Nat(1)
+	}
+	// x0 > y0 > 1
+
+	// Determine maximum result length.
+	maxLen := int(x0.Log2() - y0.Log2() + 1)
+
+	// In the following, each number x is represented
+	// as a mantissa x.m and an exponent x.e such that
+	// x = xm * 2^x.e.
+	x := fpNat{x0, 0}
+	y := fpNat{y0, 0}
+
+	// Determine a scale factor f = 2^e such that
+	// 0.5 <= y/f == y*(2^-e) < 1.0
+	// and scale y accordingly.
+	e := int(y.m.Log2()) + 1
+	y.e -= e
+
+	// t1
+	var c = 2.9142
+	const n = 14
+	t1 := fpNat{Nat(uint64(c * (1 << n))), -n}
+
+	// Compute initial value r0 for the reciprocal of y/f.
+	// r0 = t1 - 2*y
+	r := t1.sub(y.mul2())
+	two := fpNat{Nat(2), 0}
+
+	// Newton-Raphson iteration
+	p := Nat(0)
+	for i := 0; ; i++ {
+		// check if we are done
+		// TODO: Need to come up with a better test here
+		//       as it will reduce computation time significantly.
+		// q = x*r/f
+		q := x.mul(r)
+		q.e -= e
+		res := q.mant()
+		if res.Cmp(p) == 0 {
+			return res
+		}
+		p = res
+
+		// r' = r*(2 - y*r)
+		r = r.mul(two.sub(y.mul(r)))
+
+		// reduce mantissa size
+		// TODO: Find smaller bound as it will reduce
+		//       computation time massively.
+		d := int(r.m.Log2()+1) - maxLen
+		if d > 0 {
+			r = fpNat{r.m.Shr(uint(d)), r.e + d}
+		}
+	}
+
+	panic("unreachable")
+	return nil
+}
+
+
+func nrdiv(x, y Natural) (q, r Natural) {
+	q = nrDivEst(x, y)
+	r = x.Sub(y.Mul(q))
+	// if r is too large, correct q and r
+	// (usually one iteration)
+	for r.Cmp(y) >= 0 {
+		q = q.Add(Nat(1))
+		r = r.Sub(y)
+	}
+	return
+}
+
+
+func div(t *testing.T, x, y Natural) {
+	q, r := nrdiv(x, y)
+	qx, rx := x.DivMod(y)
+	if q.Cmp(qx) != 0 {
+		t.Errorf("x = %s, y = %s, got q = %s, want q = %s", x, y, q, qx)
+	}
+	if r.Cmp(rx) != 0 {
+		t.Errorf("x = %s, y = %s, got r = %s, want r = %s", x, y, r, rx)
+	}
+}
+
+
+func idiv(t *testing.T, x0, y0 uint64) { div(t, Nat(x0), Nat(y0)) }
+
+
+func TestNRDiv(t *testing.T) {
+	idiv(t, 17, 18)
+	idiv(t, 17, 17)
+	idiv(t, 17, 1)
+	idiv(t, 17, 16)
+	idiv(t, 17, 10)
+	idiv(t, 17, 9)
+	idiv(t, 17, 8)
+	idiv(t, 17, 5)
+	idiv(t, 17, 3)
+	idiv(t, 1025, 512)
+	idiv(t, 7489595, 2)
+	idiv(t, 5404679459, 78495)
+	idiv(t, 7484890589595, 7484890589594)
+	div(t, Fact(100), Fact(91))
+	div(t, Fact(1000), Fact(991))
+	//div(t, Fact(10000), Fact(9991));  // takes too long - disabled for now
+}
diff --git a/src/pkg/exp/bignum/rational.go b/src/pkg/exp/bignum/rational.go
new file mode 100644
index 000000000..378585e5f
--- /dev/null
+++ b/src/pkg/exp/bignum/rational.go
@@ -0,0 +1,205 @@
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Rational numbers
+
+package bignum
+
+import "fmt"
+
+
+// Rational represents a quotient a/b of arbitrary precision.
+//
+type Rational struct {
+	a *Integer // numerator
+	b Natural  // denominator
+}
+
+
+// MakeRat makes a rational number given a numerator a and a denominator b.
+//
+func MakeRat(a *Integer, b Natural) *Rational {
+	f := a.mant.Gcd(b) // f > 0
+	if f.Cmp(Nat(1)) != 0 {
+		a = MakeInt(a.sign, a.mant.Div(f))
+		b = b.Div(f)
+	}
+	return &Rational{a, b}
+}
+
+
+// Rat creates a small rational number with value a0/b0.
+//
+func Rat(a0 int64, b0 int64) *Rational {
+	a, b := Int(a0), Int(b0)
+	if b.sign {
+		a = a.Neg()
+	}
+	return MakeRat(a, b.mant)
+}
+
+
+// Value returns the numerator and denominator of x.
+//
+func (x *Rational) Value() (numerator *Integer, denominator Natural) {
+	return x.a, x.b
+}
+
+
+// Predicates
+
+// IsZero returns true iff x == 0.
+//
+func (x *Rational) IsZero() bool { return x.a.IsZero() }
+
+
+// IsNeg returns true iff x < 0.
+//
+func (x *Rational) IsNeg() bool { return x.a.IsNeg() }
+
+
+// IsPos returns true iff x > 0.
+//
+func (x *Rational) IsPos() bool { return x.a.IsPos() }
+
+
+// IsInt returns true iff x can be written with a denominator 1
+// in the form x == x'/1; i.e., if x is an integer value.
+//
+func (x *Rational) IsInt() bool { return x.b.Cmp(Nat(1)) == 0 }
+
+
+// Operations
+
+// Neg returns the negated value of x.
+//
+func (x *Rational) Neg() *Rational { return MakeRat(x.a.Neg(), x.b) }
+
+
+// Add returns the sum x + y.
+//
+func (x *Rational) Add(y *Rational) *Rational {
+	return MakeRat((x.a.MulNat(y.b)).Add(y.a.MulNat(x.b)), x.b.Mul(y.b))
+}
+
+
+// Sub returns the difference x - y.
+//
+func (x *Rational) Sub(y *Rational) *Rational {
+	return MakeRat((x.a.MulNat(y.b)).Sub(y.a.MulNat(x.b)), x.b.Mul(y.b))
+}
+
+
+// Mul returns the product x * y.
+//
+func (x *Rational) Mul(y *Rational) *Rational { return MakeRat(x.a.Mul(y.a), x.b.Mul(y.b)) }
+
+
+// Quo returns the quotient x / y for y != 0.
+// If y == 0, a division-by-zero run-time error occurs.
+//
+func (x *Rational) Quo(y *Rational) *Rational {
+	a := x.a.MulNat(y.b)
+	b := y.a.MulNat(x.b)
+	if b.IsNeg() {
+		a = a.Neg()
+	}
+	return MakeRat(a, b.mant)
+}
+
+
+// Cmp compares x and y. The result is an int value
+//
+//   <  0 if x <  y
+//   == 0 if x == y
+//   >  0 if x >  y
+//
+func (x *Rational) Cmp(y *Rational) int { return (x.a.MulNat(y.b)).Cmp(y.a.MulNat(x.b)) }
+
+
+// ToString converts x to a string for a given base, with 2 <= base <= 16.
+// The string representation is of the form "n" if x is an integer; otherwise
+// it is of form "n/d".
+//
+func (x *Rational) ToString(base uint) string {
+	s := x.a.ToString(base)
+	if !x.IsInt() {
+		s += "/" + x.b.ToString(base)
+	}
+	return s
+}
+
+
+// String converts x to its decimal string representation.
+// x.String() is the same as x.ToString(10).
+//
+func (x *Rational) String() string { return x.ToString(10) }
+
+
+// Format is a support routine for fmt.Formatter. It accepts
+// the formats 'b' (binary), 'o' (octal), and 'x' (hexadecimal).
+//
+func (x *Rational) Format(h fmt.State, c int) { fmt.Fprintf(h, "%s", x.ToString(fmtbase(c))) }
+
+
+// RatFromString returns the rational number corresponding to the
+// longest possible prefix of s representing a rational number in a
+// given conversion base, the actual conversion base used, and the
+// prefix length. The syntax of a rational number is:
+//
+//	rational = mantissa [exponent] .
+//	mantissa = integer ('/' natural | '.' natural) .
+//	exponent = ('e'|'E') integer .
+//
+// If the base argument is 0, the string prefix determines the actual
+// conversion base for the mantissa. A prefix of ``0x'' or ``0X'' selects
+// base 16; the ``0'' prefix selects base 8. Otherwise the selected base is 10.
+// If the mantissa is represented via a division, both the numerator and
+// denominator may have different base prefixes; in that case the base of
+// of the numerator is returned. If the mantissa contains a decimal point,
+// the base for the fractional part is the same as for the part before the
+// decimal point and the fractional part does not accept a base prefix.
+// The base for the exponent is always 10.
+//
+func RatFromString(s string, base uint) (*Rational, uint, int) {
+	// read numerator
+	a, abase, alen := IntFromString(s, base)
+	b := Nat(1)
+
+	// read denominator or fraction, if any
+	var blen int
+	if alen < len(s) {
+		ch := s[alen]
+		if ch == '/' {
+			alen++
+			b, base, blen = NatFromString(s[alen:], base)
+		} else if ch == '.' {
+			alen++
+			b, base, blen = NatFromString(s[alen:], abase)
+			assert(base == abase)
+			f := Nat(uint64(base)).Pow(uint(blen))
+			a = MakeInt(a.sign, a.mant.Mul(f).Add(b))
+			b = f
+		}
+	}
+
+	// read exponent, if any
+	rlen := alen + blen
+	if rlen < len(s) {
+		ch := s[rlen]
+		if ch == 'e' || ch == 'E' {
+			rlen++
+			e, _, elen := IntFromString(s[rlen:], 10)
+			rlen += elen
+			m := Nat(10).Pow(uint(e.mant.Value()))
+			if e.sign {
+				b = b.Mul(m)
+			} else {
+				a = a.MulNat(m)
+			}
+		}
+	}
+
+	return MakeRat(a, b), base, rlen
+}
diff --git a/src/pkg/exp/eval/eval_test.go b/src/pkg/exp/eval/eval_test.go
index 837c4fabd..1dfdfe1fd 100644
--- a/src/pkg/exp/eval/eval_test.go
+++ b/src/pkg/exp/eval/eval_test.go
@@ -5,7 +5,7 @@
 package eval
 
 import (
-	"bignum"
+	"exp/bignum"
 	"flag"
 	"fmt"
 	"log"
diff --git a/src/pkg/exp/eval/expr.go b/src/pkg/exp/eval/expr.go
index 81e9ffa93..ea8117d06 100644
--- a/src/pkg/exp/eval/expr.go
+++ b/src/pkg/exp/eval/expr.go
@@ -5,7 +5,7 @@
 package eval
 
 import (
-	"bignum"
+	"exp/bignum"
 	"fmt"
 	"go/ast"
 	"go/token"
diff --git a/src/pkg/exp/eval/expr1.go b/src/pkg/exp/eval/expr1.go
index 0e83053f4..f0a78ac4d 100644
--- a/src/pkg/exp/eval/expr1.go
+++ b/src/pkg/exp/eval/expr1.go
@@ -4,7 +4,7 @@
 package eval
 
 import (
-	"bignum"
+	"exp/bignum"
 	"log"
 )
 
diff --git a/src/pkg/exp/eval/expr_test.go b/src/pkg/exp/eval/expr_test.go
index 12914fbd5..7efa2069d 100644
--- a/src/pkg/exp/eval/expr_test.go
+++ b/src/pkg/exp/eval/expr_test.go
@@ -5,7 +5,7 @@
 package eval
 
 import (
-	"bignum"
+	"exp/bignum"
 	"testing"
 )
 
diff --git a/src/pkg/exp/eval/stmt.go b/src/pkg/exp/eval/stmt.go
index bb080375a..bcd81f04c 100644
--- a/src/pkg/exp/eval/stmt.go
+++ b/src/pkg/exp/eval/stmt.go
@@ -5,7 +5,7 @@
 package eval
 
 import (
-	"bignum"
+	"exp/bignum"
 	"log"
 	"go/ast"
 	"go/token"
diff --git a/src/pkg/exp/eval/type.go b/src/pkg/exp/eval/type.go
index 8a0a2cf2f..b0fbe2156 100644
--- a/src/pkg/exp/eval/type.go
+++ b/src/pkg/exp/eval/type.go
@@ -5,7 +5,7 @@
 package eval
 
 import (
-	"bignum"
+	"exp/bignum"
 	"go/ast"
 	"go/token"
 	"log"
diff --git a/src/pkg/exp/eval/util.go b/src/pkg/exp/eval/util.go
index 6508346dd..ffe13e170 100644
--- a/src/pkg/exp/eval/util.go
+++ b/src/pkg/exp/eval/util.go
@@ -5,7 +5,7 @@
 package eval
 
 import (
-	"bignum"
+	"exp/bignum"
 )
 
 // TODO(austin): Maybe add to bignum in more general form
diff --git a/src/pkg/exp/eval/value.go b/src/pkg/exp/eval/value.go
index 153349c43..dce4bfcf3 100644
--- a/src/pkg/exp/eval/value.go
+++ b/src/pkg/exp/eval/value.go
@@ -5,7 +5,7 @@
 package eval
 
 import (
-	"bignum"
+	"exp/bignum"
 	"fmt"
 )