libgo: Update to current sources.

git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@192704 138bc75d-0d04-0410-961f-82ee72b054a4

23 files changed, 815 insertions, 320 deletions

diff --git a/libgo/go/math/acosh.go b/libgo/go/math/acosh.go
index c6c8645e1ae..e394008b078 100644
--- a/libgo/go/math/acosh.go
+++ b/libgo/go/math/acosh.go
@@ -33,7 +33,7 @@ package math
 //	acosh(NaN) is NaN without signal.
 //
 
-// Acosh(x) calculates the inverse hyperbolic cosine of x.
+// Acosh returns the inverse hyperbolic cosine of x.
 //
 // Special cases are:
 //	Acosh(+Inf) = +Inf
diff --git a/libgo/go/math/all_test.go b/libgo/go/math/all_test.go
index 35c33ce3859..cdea8035f9a 100644
--- a/libgo/go/math/all_test.go
+++ b/libgo/go/math/all_test.go
@@ -1128,11 +1128,11 @@ var vfgammaSC = []float64{
 	NaN(),
 }
 var gammaSC = []float64{
+	NaN(),
+	NaN(),
 	Inf(-1),
 	Inf(1),
 	Inf(1),
-	Inf(1),
-	Inf(1),
 	NaN(),
 }
 
diff --git a/libgo/go/math/asinh.go b/libgo/go/math/asinh.go
index 0defbb9bef2..ff2de0215f5 100644
--- a/libgo/go/math/asinh.go
+++ b/libgo/go/math/asinh.go
@@ -30,7 +30,7 @@ package math
 //	         := sign(x)*log1p(|x| + x**2/(1 + sqrt(1+x**2)))
 //
 
-// Asinh(x) calculates the inverse hyperbolic sine of x.
+// Asinh returns the inverse hyperbolic sine of x.
 //
 // Special cases are:
 //	Asinh(±0) = ±0
diff --git a/libgo/go/math/atan.go b/libgo/go/math/atan.go
index b739721e81c..986ba068661 100644
--- a/libgo/go/math/atan.go
+++ b/libgo/go/math/atan.go
@@ -6,44 +6,85 @@ package math
 
 /*
 	Floating-point arctangent.
-
-	Atan returns the value of the arctangent of its
-	argument in the range [-pi/2,pi/2].
-	There are no error returns.
-	Coefficients are #5077 from Hart & Cheney. (19.56D)
 */
 
-// xatan evaluates a series valid in the
-// range [-0.414...,+0.414...]. (tan(pi/8))
-func xatan(arg float64) float64 {
+// The original C code, the long comment, and the constants below were
+// from http://netlib.sandia.gov/cephes/cmath/atan.c, available from
+// http://www.netlib.org/cephes/cmath.tgz.
+// The go code is a version of the original C.
+//
+// atan.c
+// Inverse circular tangent (arctangent)
+//
+// SYNOPSIS:
+// double x, y, atan();
+// y = atan( x );
+//
+// DESCRIPTION:
+// Returns radian angle between -pi/2 and +pi/2 whose tangent is x.
+//
+// Range reduction is from three intervals into the interval from zero to 0.66.
+// The approximant uses a rational function of degree 4/5 of the form
+// x + x**3 P(x)/Q(x).
+//
+// ACCURACY:
+//                      Relative error:
+// arithmetic   domain    # trials  peak     rms
+//    DEC       -10, 10   50000     2.4e-17  8.3e-18
+//    IEEE      -10, 10   10^6      1.8e-16  5.0e-17
+//
+// Cephes Math Library Release 2.8:  June, 2000
+// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+//
+// The readme file at http://netlib.sandia.gov/cephes/ says:
+//    Some software in this archive may be from the book _Methods and
+// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
+// International, 1989) or from the Cephes Mathematical Library, a
+// commercial product. In either event, it is copyrighted by the author.
+// What you see here may be used freely but it comes with no support or
+// guarantee.
+//
+//   The two known misprints in the book are repaired here in the
+// source listings for the gamma function and the incomplete beta
+// integral.
+//
+//   Stephen L. Moshier
+//   moshier@na-net.ornl.gov
+
+// xatan evaluates a series valid in the range [0, 0.66].
+func xatan(x float64) float64 {
 	const (
-		P4 = .161536412982230228262e2
-		P3 = .26842548195503973794141e3
-		P2 = .11530293515404850115428136e4
-		P1 = .178040631643319697105464587e4
-		P0 = .89678597403663861959987488e3
-		Q4 = .5895697050844462222791e2
-		Q3 = .536265374031215315104235e3
-		Q2 = .16667838148816337184521798e4
-		Q1 = .207933497444540981287275926e4
-		Q0 = .89678597403663861962481162e3
+		P0 = -8.750608600031904122785e-01
+		P1 = -1.615753718733365076637e+01
+		P2 = -7.500855792314704667340e+01
+		P3 = -1.228866684490136173410e+02
+		P4 = -6.485021904942025371773e+01
+		Q0 = +2.485846490142306297962e+01
+		Q1 = +1.650270098316988542046e+02
+		Q2 = +4.328810604912902668951e+02
+		Q3 = +4.853903996359136964868e+02
+		Q4 = +1.945506571482613964425e+02
 	)
-	sq := arg * arg
-	value := ((((P4*sq+P3)*sq+P2)*sq+P1)*sq + P0)
-	value = value / (((((sq+Q4)*sq+Q3)*sq+Q2)*sq+Q1)*sq + Q0)
-	return value * arg
+	z := x * x
+	z = z * ((((P0*z+P1)*z+P2)*z+P3)*z + P4) / (((((z+Q0)*z+Q1)*z+Q2)*z+Q3)*z + Q4)
+	z = x*z + x
+	return z
 }
 
 // satan reduces its argument (known to be positive)
-// to the range [0,0.414...] and calls xatan.
-func satan(arg float64) float64 {
-	if arg < Sqrt2-1 {
-		return xatan(arg)
+// to the range [0, 0.66] and calls xatan.
+func satan(x float64) float64 {
+	const (
+		Morebits = 6.123233995736765886130e-17 // pi/2 = PIO2 + Morebits
+		Tan3pio8 = 2.41421356237309504880      // tan(3*pi/8)
+	)
+	if x <= 0.66 {
+		return xatan(x)
 	}
-	if arg > Sqrt2+1 {
-		return Pi/2 - xatan(1/arg)
+	if x > Tan3pio8 {
+		return Pi/2 - xatan(1/x) + Morebits
 	}
-	return Pi/4 + xatan((arg-1)/(arg+1))
+	return Pi/4 + xatan((x-1)/(x+1)) + 0.5*Morebits
 }
 
 // Atan returns the arctangent of x.
diff --git a/libgo/go/math/atanh.go b/libgo/go/math/atanh.go
index 5b5d468559b..113d5c103c0 100644
--- a/libgo/go/math/atanh.go
+++ b/libgo/go/math/atanh.go
@@ -36,7 +36,7 @@ package math
 //	atanh(+-1) is +-INF with signal.
 //
 
-// Atanh(x) calculates the inverse hyperbolic tangent of x.
+// Atanh returns the inverse hyperbolic tangent of x.
 //
 // Special cases are:
 //	Atanh(1) = +Inf
diff --git a/libgo/go/math/big/arith_test.go b/libgo/go/math/big/arith_test.go
index c7e3d284c2d..3615a659c3b 100644
--- a/libgo/go/math/big/arith_test.go
+++ b/libgo/go/math/big/arith_test.go
@@ -4,7 +4,10 @@
 
 package big
 
-import "testing"
+import (
+	"math/rand"
+	"testing"
+)
 
 type funWW func(x, y, c Word) (z1, z0 Word)
 type argWW struct {
@@ -100,6 +103,43 @@ func TestFunVV(t *testing.T) {
 	}
 }
 
+// Always the same seed for reproducible results.
+var rnd = rand.New(rand.NewSource(0))
+
+func rndW() Word {
+	return Word(rnd.Int63()<<1 | rnd.Int63n(2))
+}
+
+func rndV(n int) []Word {
+	v := make([]Word, n)
+	for i := range v {
+		v[i] = rndW()
+	}
+	return v
+}
+
+func benchmarkFunVV(b *testing.B, f funVV, n int) {
+	x := rndV(n)
+	y := rndV(n)
+	z := make([]Word, n)
+	b.SetBytes(int64(n * _W))
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		f(z, x, y)
+	}
+}
+
+func BenchmarkAddVV_1(b *testing.B)   { benchmarkFunVV(b, addVV, 1) }
+func BenchmarkAddVV_2(b *testing.B)   { benchmarkFunVV(b, addVV, 2) }
+func BenchmarkAddVV_3(b *testing.B)   { benchmarkFunVV(b, addVV, 3) }
+func BenchmarkAddVV_4(b *testing.B)   { benchmarkFunVV(b, addVV, 4) }
+func BenchmarkAddVV_5(b *testing.B)   { benchmarkFunVV(b, addVV, 5) }
+func BenchmarkAddVV_1e1(b *testing.B) { benchmarkFunVV(b, addVV, 1e1) }
+func BenchmarkAddVV_1e2(b *testing.B) { benchmarkFunVV(b, addVV, 1e2) }
+func BenchmarkAddVV_1e3(b *testing.B) { benchmarkFunVV(b, addVV, 1e3) }
+func BenchmarkAddVV_1e4(b *testing.B) { benchmarkFunVV(b, addVV, 1e4) }
+func BenchmarkAddVV_1e5(b *testing.B) { benchmarkFunVV(b, addVV, 1e5) }
+
 type funVW func(z, x []Word, y Word) (c Word)
 type argVW struct {
 	z, x nat
@@ -210,6 +250,28 @@ func TestFunVW(t *testing.T) {
 	}
 }
 
+func benchmarkFunVW(b *testing.B, f funVW, n int) {
+	x := rndV(n)
+	y := rndW()
+	z := make([]Word, n)
+	b.SetBytes(int64(n * _W))
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		f(z, x, y)
+	}
+}
+
+func BenchmarkAddVW_1(b *testing.B)   { benchmarkFunVW(b, addVW, 1) }
+func BenchmarkAddVW_2(b *testing.B)   { benchmarkFunVW(b, addVW, 2) }
+func BenchmarkAddVW_3(b *testing.B)   { benchmarkFunVW(b, addVW, 3) }
+func BenchmarkAddVW_4(b *testing.B)   { benchmarkFunVW(b, addVW, 4) }
+func BenchmarkAddVW_5(b *testing.B)   { benchmarkFunVW(b, addVW, 5) }
+func BenchmarkAddVW_1e1(b *testing.B) { benchmarkFunVW(b, addVW, 1e1) }
+func BenchmarkAddVW_1e2(b *testing.B) { benchmarkFunVW(b, addVW, 1e2) }
+func BenchmarkAddVW_1e3(b *testing.B) { benchmarkFunVW(b, addVW, 1e3) }
+func BenchmarkAddVW_1e4(b *testing.B) { benchmarkFunVW(b, addVW, 1e4) }
+func BenchmarkAddVW_1e5(b *testing.B) { benchmarkFunVW(b, addVW, 1e5) }
+
 type funVWW func(z, x []Word, y, r Word) (c Word)
 type argVWW struct {
 	z, x nat
@@ -334,6 +396,28 @@ func TestMulAddWWW(t *testing.T) {
 	}
 }
 
+func benchmarkAddMulVVW(b *testing.B, n int) {
+	x := rndV(n)
+	y := rndW()
+	z := make([]Word, n)
+	b.SetBytes(int64(n * _W))
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		addMulVVW(z, x, y)
+	}
+}
+
+func BenchmarkAddMulVVW_1(b *testing.B)   { benchmarkAddMulVVW(b, 1) }
+func BenchmarkAddMulVVW_2(b *testing.B)   { benchmarkAddMulVVW(b, 2) }
+func BenchmarkAddMulVVW_3(b *testing.B)   { benchmarkAddMulVVW(b, 3) }
+func BenchmarkAddMulVVW_4(b *testing.B)   { benchmarkAddMulVVW(b, 4) }
+func BenchmarkAddMulVVW_5(b *testing.B)   { benchmarkAddMulVVW(b, 5) }
+func BenchmarkAddMulVVW_1e1(b *testing.B) { benchmarkAddMulVVW(b, 1e1) }
+func BenchmarkAddMulVVW_1e2(b *testing.B) { benchmarkAddMulVVW(b, 1e2) }
+func BenchmarkAddMulVVW_1e3(b *testing.B) { benchmarkAddMulVVW(b, 1e3) }
+func BenchmarkAddMulVVW_1e4(b *testing.B) { benchmarkAddMulVVW(b, 1e4) }
+func BenchmarkAddMulVVW_1e5(b *testing.B) { benchmarkAddMulVVW(b, 1e5) }
+
 func testWordBitLen(t *testing.T, fname string, f func(Word) int) {
 	for i := 0; i <= _W; i++ {
 		x := Word(1) << uint(i-1) // i == 0 => x == 0
diff --git a/libgo/go/math/big/calibrate_test.go b/libgo/go/math/big/calibrate_test.go
index efe1837bba3..f69ffbf5cf1 100644
--- a/libgo/go/math/big/calibrate_test.go
+++ b/libgo/go/math/big/calibrate_test.go
@@ -21,15 +21,17 @@ import (
 
 var calibrate = flag.Bool("calibrate", false, "run calibration test")
 
-// measure returns the time to run f
-func measure(f func()) time.Duration {
-	const N = 100
-	start := time.Now()
-	for i := N; i > 0; i-- {
-		f()
-	}
-	stop := time.Now()
-	return stop.Sub(start) / N
+func karatsubaLoad(b *testing.B) {
+	BenchmarkMul(b)
+}
+
+// measureKaratsuba returns the time to run a Karatsuba-relevant benchmark
+// given Karatsuba threshold th.
+func measureKaratsuba(th int) time.Duration {
+	th, karatsubaThreshold = karatsubaThreshold, th
+	res := testing.Benchmark(karatsubaLoad)
+	karatsubaThreshold = th
+	return time.Duration(res.NsPerOp())
 }
 
 func computeThresholds() {
@@ -37,35 +39,33 @@ func computeThresholds() {
 	fmt.Printf("(run repeatedly for good results)\n")
 
 	// determine Tk, the work load execution time using basic multiplication
-	karatsubaThreshold = 1e9 // disable karatsuba
-	Tb := measure(benchmarkMulLoad)
-	fmt.Printf("Tb = %dns\n", Tb)
+	Tb := measureKaratsuba(1e9) // th == 1e9 => Karatsuba multiplication disabled
+	fmt.Printf("Tb = %10s\n", Tb)
 
 	// thresholds
-	n := 8 // any lower values for the threshold lead to very slow multiplies
+	th := 4
 	th1 := -1
 	th2 := -1
 
 	var deltaOld time.Duration
-	for count := -1; count != 0; count-- {
+	for count := -1; count != 0 && th < 128; count-- {
 		// determine Tk, the work load execution time using Karatsuba multiplication
-		karatsubaThreshold = n // enable karatsuba
-		Tk := measure(benchmarkMulLoad)
+		Tk := measureKaratsuba(th)
 
 		// improvement over Tb
 		delta := (Tb - Tk) * 100 / Tb
 
-		fmt.Printf("n = %3d  Tk = %8dns  %4d%%", n, Tk, delta)
+		fmt.Printf("th = %3d  Tk = %10s  %4d%%", th, Tk, delta)
 
 		// determine break-even point
 		if Tk < Tb && th1 < 0 {
-			th1 = n
+			th1 = th
 			fmt.Print("  break-even point")
 		}
 
 		// determine diminishing return
 		if 0 < delta && delta < deltaOld && th2 < 0 {
-			th2 = n
+			th2 = th
 			fmt.Print("  diminishing return")
 		}
 		deltaOld = delta
@@ -74,10 +74,10 @@ func computeThresholds() {
 
 		// trigger counter
 		if th1 >= 0 && th2 >= 0 && count < 0 {
-			count = 20 // this many extra measurements after we got both thresholds
+			count = 10 // this many extra measurements after we got both thresholds
 		}
 
-		n++
+		th++
 	}
 }
 
diff --git a/libgo/go/math/big/gcd_test.go b/libgo/go/math/big/gcd_test.go
new file mode 100644
index 00000000000..c0b9f583000
--- /dev/null
+++ b/libgo/go/math/big/gcd_test.go
@@ -0,0 +1,47 @@
+// Copyright 2012 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 a GCD benchmark.
+// Usage: go test math/big -test.bench GCD
+
+package big
+
+import (
+	"math/rand"
+	"testing"
+)
+
+// randInt returns a pseudo-random Int in the range [1<<(size-1), (1<<size) - 1]
+func randInt(r *rand.Rand, size uint) *Int {
+	n := new(Int).Lsh(intOne, size-1)
+	x := new(Int).Rand(r, n)
+	return x.Add(x, n) // make sure result > 1<<(size-1)
+}
+
+func runGCD(b *testing.B, aSize, bSize uint) {
+	b.StopTimer()
+	var r = rand.New(rand.NewSource(1234))
+	aa := randInt(r, aSize)
+	bb := randInt(r, bSize)
+	b.StartTimer()
+	for i := 0; i < b.N; i++ {
+		new(Int).GCD(nil, nil, aa, bb)
+	}
+}
+
+func BenchmarkGCD10x10(b *testing.B)         { runGCD(b, 10, 10) }
+func BenchmarkGCD10x100(b *testing.B)        { runGCD(b, 10, 100) }
+func BenchmarkGCD10x1000(b *testing.B)       { runGCD(b, 10, 1000) }
+func BenchmarkGCD10x10000(b *testing.B)      { runGCD(b, 10, 10000) }
+func BenchmarkGCD10x100000(b *testing.B)     { runGCD(b, 10, 100000) }
+func BenchmarkGCD100x100(b *testing.B)       { runGCD(b, 100, 100) }
+func BenchmarkGCD100x1000(b *testing.B)      { runGCD(b, 100, 1000) }
+func BenchmarkGCD100x10000(b *testing.B)     { runGCD(b, 100, 10000) }
+func BenchmarkGCD100x100000(b *testing.B)    { runGCD(b, 100, 100000) }
+func BenchmarkGCD1000x1000(b *testing.B)     { runGCD(b, 1000, 1000) }
+func BenchmarkGCD1000x10000(b *testing.B)    { runGCD(b, 1000, 10000) }
+func BenchmarkGCD1000x100000(b *testing.B)   { runGCD(b, 1000, 100000) }
+func BenchmarkGCD10000x10000(b *testing.B)   { runGCD(b, 10000, 10000) }
+func BenchmarkGCD10000x100000(b *testing.B)  { runGCD(b, 10000, 100000) }
+func BenchmarkGCD100000x100000(b *testing.B) { runGCD(b, 100000, 100000) }
diff --git a/libgo/go/math/big/int.go b/libgo/go/math/big/int.go
index cd2cd0e2da0..4bd7958ae50 100644
--- a/libgo/go/math/big/int.go
+++ b/libgo/go/math/big/int.go
@@ -581,12 +581,12 @@ func (z *Int) Exp(x, y, m *Int) *Int {
 	return z
 }
 
-// GCD sets z to the greatest common divisor of a and b, which must be
-// positive numbers, and returns z.
+// GCD sets z to the greatest common divisor of a and b, which both must
+// be > 0, and returns z.
 // If x and y are not nil, GCD sets x and y such that z = a*x + b*y.
-// If either a or b is not positive, GCD sets z = x = y = 0.
+// If either a or b is <= 0, GCD sets z = x = y = 0.
 func (z *Int) GCD(x, y, a, b *Int) *Int {
-	if a.neg || b.neg {
+	if a.Sign() <= 0 || b.Sign() <= 0 {
 		z.SetInt64(0)
 		if x != nil {
 			x.SetInt64(0)
@@ -596,6 +596,9 @@ func (z *Int) GCD(x, y, a, b *Int) *Int {
 		}
 		return z
 	}
+	if x == nil && y == nil {
+		return z.binaryGCD(a, b)
+	}
 
 	A := new(Int).Set(a)
 	B := new(Int).Set(b)
@@ -640,6 +643,63 @@ func (z *Int) GCD(x, y, a, b *Int) *Int {
 	return z
 }
 
+// binaryGCD sets z to the greatest common divisor of a and b, which both must
+// be > 0, and returns z.
+// See Knuth, The Art of Computer Programming, Vol. 2, Section 4.5.2, Algorithm B.
+func (z *Int) binaryGCD(a, b *Int) *Int {
+	u := z
+	v := new(Int)
+
+	// use one Euclidean iteration to ensure that u and v are approx. the same size
+	switch {
+	case len(a.abs) > len(b.abs):
+		u.Set(b)
+		v.Rem(a, b)
+	case len(a.abs) < len(b.abs):
+		u.Set(a)
+		v.Rem(b, a)
+	default:
+		u.Set(a)
+		v.Set(b)
+	}
+
+	// v might be 0 now
+	if len(v.abs) == 0 {
+		return u
+	}
+	// u > 0 && v > 0
+
+	// determine largest k such that u = u' << k, v = v' << k
+	k := u.abs.trailingZeroBits()
+	if vk := v.abs.trailingZeroBits(); vk < k {
+		k = vk
+	}
+	u.Rsh(u, k)
+	v.Rsh(v, k)
+
+	// determine t (we know that u > 0)
+	t := new(Int)
+	if u.abs[0]&1 != 0 {
+		// u is odd
+		t.Neg(v)
+	} else {
+		t.Set(u)
+	}
+
+	for len(t.abs) > 0 {
+		// reduce t
+		t.Rsh(t, t.abs.trailingZeroBits())
+		if t.neg {
+			v.Neg(t)
+		} else {
+			u.Set(t)
+		}
+		t.Sub(u, v)
+	}
+
+	return u.Lsh(u, k)
+}
+
 // ProbablyPrime performs n Miller-Rabin tests to check whether x is prime.
 // If it returns true, x is prime with probability 1 - 1/4^n.
 // If it returns false, x is not prime.
@@ -697,6 +757,13 @@ func (z *Int) Rsh(x *Int, n uint) *Int {
 // Bit returns the value of the i'th bit of x. That is, it
 // returns (x>>i)&1. The bit index i must be >= 0.
 func (x *Int) Bit(i int) uint {
+	if i == 0 {
+		// optimization for common case: odd/even test of x
+		if len(x.abs) > 0 {
+			return uint(x.abs[0] & 1) // bit 0 is same for -x
+		}
+		return 0
+	}
 	if i < 0 {
 		panic("negative bit index")
 	}
@@ -894,3 +961,19 @@ func (z *Int) GobDecode(buf []byte) error {
 	z.abs = z.abs.setBytes(buf[1:])
 	return nil
 }
+
+// MarshalJSON implements the json.Marshaler interface.
+func (x *Int) MarshalJSON() ([]byte, error) {
+	// TODO(gri): get rid of the []byte/string conversions
+	return []byte(x.String()), nil
+}
+
+// UnmarshalJSON implements the json.Unmarshaler interface.
+func (z *Int) UnmarshalJSON(x []byte) error {
+	// TODO(gri): get rid of the []byte/string conversions
+	_, ok := z.SetString(string(x), 0)
+	if !ok {
+		return fmt.Errorf("math/big: cannot unmarshal %s into a *big.Int", x)
+	}
+	return nil
+}
diff --git a/libgo/go/math/big/int_test.go b/libgo/go/math/big/int_test.go
index 9700a9b5a7c..27834cec6af 100644
--- a/libgo/go/math/big/int_test.go
+++ b/libgo/go/math/big/int_test.go
@@ -8,6 +8,7 @@ import (
 	"bytes"
 	"encoding/gob"
 	"encoding/hex"
+	"encoding/json"
 	"fmt"
 	"math/rand"
 	"testing"
@@ -817,14 +818,12 @@ func TestExp(t *testing.T) {
 }
 
 func checkGcd(aBytes, bBytes []byte) bool {
-	a := new(Int).SetBytes(aBytes)
-	b := new(Int).SetBytes(bBytes)
-
 	x := new(Int)
 	y := new(Int)
-	d := new(Int)
+	a := new(Int).SetBytes(aBytes)
+	b := new(Int).SetBytes(bBytes)
 
-	d.GCD(x, y, a, b)
+	d := new(Int).GCD(x, y, a, b)
 	x.Mul(x, a)
 	y.Mul(y, b)
 	x.Add(x, y)
@@ -833,32 +832,70 @@ func checkGcd(aBytes, bBytes []byte) bool {
 }
 
 var gcdTests = []struct {
-	a, b    int64
-	d, x, y int64
+	d, x, y, a, b string
 }{
-	{120, 23, 1, -9, 47},
-}
-
-func TestGcd(t *testing.T) {
-	for i, test := range gcdTests {
-		a := NewInt(test.a)
-		b := NewInt(test.b)
+	// a <= 0 || b <= 0
+	{"0", "0", "0", "0", "0"},
+	{"0", "0", "0", "0", "7"},
+	{"0", "0", "0", "11", "0"},
+	{"0", "0", "0", "-77", "35"},
+	{"0", "0", "0", "64515", "-24310"},
+	{"0", "0", "0", "-64515", "-24310"},
+
+	{"1", "-9", "47", "120", "23"},
+	{"7", "1", "-2", "77", "35"},
+	{"935", "-3", "8", "64515", "24310"},
+	{"935000000000000000", "-3", "8", "64515000000000000000", "24310000000000000000"},
+	{"1", "-221", "22059940471369027483332068679400581064239780177629666810348940098015901108344", "98920366548084643601728869055592650835572950932266967461790948584315647051443", "991"},
+
+	// test early exit (after one Euclidean iteration) in binaryGCD
+	{"1", "", "", "1", "98920366548084643601728869055592650835572950932266967461790948584315647051443"},
+}
+
+func testGcd(t *testing.T, d, x, y, a, b *Int) {
+	var X *Int
+	if x != nil {
+		X = new(Int)
+	}
+	var Y *Int
+	if y != nil {
+		Y = new(Int)
+	}
 
-		x := new(Int)
-		y := new(Int)
-		d := new(Int)
+	D := new(Int).GCD(X, Y, a, b)
+	if D.Cmp(d) != 0 {
+		t.Errorf("GCD(%s, %s): got d = %s, want %s", a, b, D, d)
+	}
+	if x != nil && X.Cmp(x) != 0 {
+		t.Errorf("GCD(%s, %s): got x = %s, want %s", a, b, X, x)
+	}
+	if y != nil && Y.Cmp(y) != 0 {
+		t.Errorf("GCD(%s, %s): got y = %s, want %s", a, b, Y, y)
+	}
 
-		expectedX := NewInt(test.x)
-		expectedY := NewInt(test.y)
-		expectedD := NewInt(test.d)
+	// binaryGCD requires a > 0 && b > 0
+	if a.Sign() <= 0 || b.Sign() <= 0 {
+		return
+	}
 
-		d.GCD(x, y, a, b)
+	D.binaryGCD(a, b)
+	if D.Cmp(d) != 0 {
+		t.Errorf("binaryGcd(%s, %s): got d = %s, want %s", a, b, D, d)
+	}
+}
 
-		if expectedX.Cmp(x) != 0 ||
-			expectedY.Cmp(y) != 0 ||
-			expectedD.Cmp(d) != 0 {
-			t.Errorf("#%d got (%s %s %s) want (%s %s %s)", i, x, y, d, expectedX, expectedY, expectedD)
-		}
+func TestGcd(t *testing.T) {
+	for _, test := range gcdTests {
+		d, _ := new(Int).SetString(test.d, 0)
+		x, _ := new(Int).SetString(test.x, 0)
+		y, _ := new(Int).SetString(test.y, 0)
+		a, _ := new(Int).SetString(test.a, 0)
+		b, _ := new(Int).SetString(test.b, 0)
+
+		testGcd(t, d, nil, nil, a, b)
+		testGcd(t, d, x, nil, a, b)
+		testGcd(t, d, nil, y, a, b)
+		testGcd(t, d, x, y, a, b)
 	}
 
 	quick.Check(checkGcd, nil)
@@ -1368,8 +1405,12 @@ func TestModInverse(t *testing.T) {
 	}
 }
 
-// used by TestIntGobEncoding and TestRatGobEncoding
-var gobEncodingTests = []string{
+var encodingTests = []string{
+	"-539345864568634858364538753846587364875430589374589",
+	"-678645873",
+	"-100",
+	"-2",
+	"-1",
 	"0",
 	"1",
 	"2",
@@ -1383,26 +1424,37 @@ func TestIntGobEncoding(t *testing.T) {
 	var medium bytes.Buffer
 	enc := gob.NewEncoder(&medium)
 	dec := gob.NewDecoder(&medium)
-	for i, test := range gobEncodingTests {
-		for j := 0; j < 2; j++ {
-			medium.Reset() // empty buffer for each test case (in case of failures)
-			stest := test
-			if j != 0 {
-				// negative numbers
-				stest = "-" + test
-			}
-			var tx Int
-			tx.SetString(stest, 10)
-			if err := enc.Encode(&tx); err != nil {
-				t.Errorf("#%d%c: encoding failed: %s", i, 'a'+j, err)
-			}
-			var rx Int
-			if err := dec.Decode(&rx); err != nil {
-				t.Errorf("#%d%c: decoding failed: %s", i, 'a'+j, err)
-			}
-			if rx.Cmp(&tx) != 0 {
-				t.Errorf("#%d%c: transmission failed: got %s want %s", i, 'a'+j, &rx, &tx)
-			}
+	for _, test := range encodingTests {
+		medium.Reset() // empty buffer for each test case (in case of failures)
+		var tx Int
+		tx.SetString(test, 10)
+		if err := enc.Encode(&tx); err != nil {
+			t.Errorf("encoding of %s failed: %s", &tx, err)
+		}
+		var rx Int
+		if err := dec.Decode(&rx); err != nil {
+			t.Errorf("decoding of %s failed: %s", &tx, err)
+		}
+		if rx.Cmp(&tx) != 0 {
+			t.Errorf("transmission of %s failed: got %s want %s", &tx, &rx, &tx)
+		}
+	}
+}
+
+func TestIntJSONEncoding(t *testing.T) {
+	for _, test := range encodingTests {
+		var tx Int
+		tx.SetString(test, 10)
+		b, err := json.Marshal(&tx)
+		if err != nil {
+			t.Errorf("marshaling of %s failed: %s", &tx, err)
+		}
+		var rx Int
+		if err := json.Unmarshal(b, &rx); err != nil {
+			t.Errorf("unmarshaling of %s failed: %s", &tx, err)
+		}
+		if rx.Cmp(&tx) != 0 {
+			t.Errorf("JSON encoding of %s failed: got %s want %s", &tx, &rx, &tx)
 		}
 	}
 }
diff --git a/libgo/go/math/big/nat.go b/libgo/go/math/big/nat.go
index 6d81823bb4a..2d5a5c9587d 100644
--- a/libgo/go/math/big/nat.go
+++ b/libgo/go/math/big/nat.go
@@ -236,7 +236,7 @@ func karatsubaSub(z, x nat, n int) {
 // Operands that are shorter than karatsubaThreshold are multiplied using
 // "grade school" multiplication; for longer operands the Karatsuba algorithm
 // is used.
-var karatsubaThreshold int = 32 // computed by calibrate.go
+var karatsubaThreshold int = 40 // computed by calibrate.go
 
 // karatsuba multiplies x and y and leaves the result in z.
 // Both x and y must have the same length n and n must be a
@@ -342,7 +342,7 @@ func alias(x, y nat) bool {
 	return cap(x) > 0 && cap(y) > 0 && &x[0:cap(x)][cap(x)-1] == &y[0:cap(y)][cap(y)-1]
 }
 
-// addAt implements z += x*(1<<(_W*i)); z must be long enough.
+// addAt implements z += x<<(_W*i); z must be long enough.
 // (we don't use nat.add because we need z to stay the same
 // slice, and we don't need to normalize z after each addition)
 func addAt(z, x nat, i int) {
@@ -405,8 +405,8 @@ func (z nat) mul(x, y nat) nat {
 
 	// determine Karatsuba length k such that
 	//
-	//   x = x1*b + x0
-	//   y = y1*b + y0  (and k <= len(y), which implies k <= len(x))
+	//   x = xh*b + x0  (0 <= x0 < b)
+	//   y = yh*b + y0  (0 <= y0 < b)
 	//   b = 1<<(_W*k)  ("base" of digits xi, yi)
 	//
 	k := karatsubaLen(n)
@@ -417,27 +417,44 @@ func (z nat) mul(x, y nat) nat {
 	y0 := y[0:k]              // y0 is not normalized
 	z = z.make(max(6*k, m+n)) // enough space for karatsuba of x0*y0 and full result of x*y
 	karatsuba(z, x0, y0)
-	z = z[0 : m+n] // z has final length but may be incomplete, upper portion is garbage
-
-	// If x1 and/or y1 are not 0, add missing terms to z explicitly:
-	//
-	//     m+n       2*k       0
-	//   z = [   ...   | x0*y0 ]
-	//     +   [ x1*y1 ]
-	//     +   [ x1*y0 ]
-	//     +   [ x0*y1 ]
+	z = z[0 : m+n]  // z has final length but may be incomplete
+	z[2*k:].clear() // upper portion of z is garbage (and 2*k <= m+n since k <= n <= m)
+
+	// If xh != 0 or yh != 0, add the missing terms to z. For
+	// 
+	//   xh = xi*b^i + ... + x2*b^2 + x1*b (0 <= xi < b) 
+	//   yh =                         y1*b (0 <= y1 < b) 
+	// 
+	// the missing terms are 
+	// 
+	//   x0*y1*b and xi*y0*b^i, xi*y1*b^(i+1) for i > 0 
+	// 
+	// since all the yi for i > 1 are 0 by choice of k: If any of them 
+	// were > 0, then yh >= b^2 and thus y >= b^2. Then k' = k*2 would 
+	// be a larger valid threshold contradicting the assumption about k. 
 	//
 	if k < n || m != n {
-		x1 := x[k:] // x1 is normalized because x is
-		y1 := y[k:] // y1 is normalized because y is
 		var t nat
-		t = t.mul(x1, y1)
-		copy(z[2*k:], t)
-		z[2*k+len(t):].clear() // upper portion of z is garbage
-		t = t.mul(x1, y0.norm())
-		addAt(z, t, k)
-		t = t.mul(x0.norm(), y1)
+
+		// add x0*y1*b
+		x0 := x0.norm()
+		y1 := y[k:]       // y1 is normalized because y is
+		t = t.mul(x0, y1) // update t so we don't lose t's underlying array
 		addAt(z, t, k)
+
+		// add xi*y0<<i, xi*y1*b<<(i+k)
+		y0 := y0.norm()
+		for i := k; i < len(x); i += k {
+			xi := x[i:]
+			if len(xi) > k {
+				xi = xi[:k]
+			}
+			xi = xi.norm()
+			t = t.mul(xi, y0)
+			addAt(z, t, i)
+			t = t.mul(xi, y1)
+			addAt(z, t, i+k)
+		}
 	}
 
 	return z.norm()
@@ -493,14 +510,9 @@ func (z nat) div(z2, u, v nat) (q, r nat) {
 	}
 
 	if len(v) == 1 {
-		var rprime Word
-		q, rprime = z.divW(u, v[0])
-		if rprime > 0 {
-			r = z2.make(1)
-			r[0] = rprime
-		} else {
-			r = z2.make(0)
-		}
+		var r2 Word
+		q, r2 = z.divW(u, v[0])
+		r = z2.setWord(r2)
 		return
 	}
 
@@ -740,7 +752,7 @@ func (x nat) string(charset string) string {
 	// convert power of two and non power of two bases separately
 	if b == b&-b {
 		// shift is base-b digit size in bits
-		shift := uint(trailingZeroBits(b)) // shift > 0 because b >= 2
+		shift := trailingZeroBits(b) // shift > 0 because b >= 2
 		mask := Word(1)<<shift - 1
 		w := x[0]
 		nbits := uint(_W) // number of unprocessed bits in w
@@ -993,10 +1005,9 @@ var deBruijn64Lookup = []byte{
 	54, 26, 40, 15, 34, 20, 31, 10, 25, 14, 19, 9, 13, 8, 7, 6,
 }
 
-// trailingZeroBits returns the number of consecutive zero bits on the right
-// side of the given Word.
-// See Knuth, volume 4, section 7.3.1
-func trailingZeroBits(x Word) int {
+// trailingZeroBits returns the number of consecutive least significant zero
+// bits of x.
+func trailingZeroBits(x Word) uint {
 	// x & -x leaves only the right-most bit set in the word. Let k be the
 	// index of that bit. Since only a single bit is set, the value is two
 	// to the power of k. Multiplying by a power of two is equivalent to
@@ -1005,18 +1016,33 @@ func trailingZeroBits(x Word) int {
 	// Therefore, if we have a left shifted version of this constant we can
 	// find by how many bits it was shifted by looking at which six bit
 	// substring ended up at the top of the word.
+	// (Knuth, volume 4, section 7.3.1)
 	switch _W {
 	case 32:
-		return int(deBruijn32Lookup[((x&-x)*deBruijn32)>>27])
+		return uint(deBruijn32Lookup[((x&-x)*deBruijn32)>>27])
 	case 64:
-		return int(deBruijn64Lookup[((x&-x)*(deBruijn64&_M))>>58])
+		return uint(deBruijn64Lookup[((x&-x)*(deBruijn64&_M))>>58])
 	default:
-		panic("Unknown word size")
+		panic("unknown word size")
 	}
 
 	return 0
 }
 
+// trailingZeroBits returns the number of consecutive least significant zero
+// bits of x.
+func (x nat) trailingZeroBits() uint {
+	if len(x) == 0 {
+		return 0
+	}
+	var i uint
+	for x[i] == 0 {
+		i++
+	}
+	// x[i] != 0
+	return i*_W + trailingZeroBits(x[i])
+}
+
 // z = x << s
 func (z nat) shl(x nat, s uint) nat {
 	m := len(x)
@@ -1169,29 +1195,6 @@ func (x nat) modW(d Word) (r Word) {
 	return divWVW(q, 0, x, d)
 }
 
-// powersOfTwoDecompose finds q and k with x = q * 1<<k and q is odd, or q and k are 0.
-func (x nat) powersOfTwoDecompose() (q nat, k int) {
-	if len(x) == 0 {
-		return x, 0
-	}
-
-	// One of the words must be non-zero by definition,
-	// so this loop will terminate with i < len(x), and
-	// i is the number of 0 words.
-	i := 0
-	for x[i] == 0 {
-		i++
-	}
-	n := trailingZeroBits(x[i]) // x[i] != 0
-
-	q = make(nat, len(x)-i)
-	shrVU(q, x[i:], uint(n))
-
-	q = q.norm()
-	k = i*_W + n
-	return
-}
-
 // random creates a random integer in [0..limit), using the space in z if
 // possible. n is the bit length of limit.
 func (z nat) random(rand *rand.Rand, limit nat, n int) nat {
@@ -1207,17 +1210,19 @@ func (z nat) random(rand *rand.Rand, limit nat, n int) nat {
 	mask := Word((1 << bitLengthOfMSW) - 1)
 
 	for {
-		for i := range z {
-			switch _W {
-			case 32:
+		switch _W {
+		case 32:
+			for i := range z {
 				z[i] = Word(rand.Uint32())
-			case 64:
+			}
+		case 64:
+			for i := range z {
 				z[i] = Word(rand.Uint32()) | Word(rand.Uint32())<<32
 			}
+		default:
+			panic("unknown word size")
 		}
-
 		z[len(limit)-1] &= mask
-
 		if z.cmp(limit) < 0 {
 			break
 		}
@@ -1230,7 +1235,7 @@ func (z nat) random(rand *rand.Rand, limit nat, n int) nat {
 // reuses the storage of z if possible.
 func (z nat) expNN(x, y, m nat) nat {
 	if alias(z, x) || alias(z, y) {
-		// We cannot allow in place modification of x or y.
+		// We cannot allow in-place modification of x or y.
 		z = nil
 	}
 
@@ -1259,15 +1264,21 @@ func (z nat) expNN(x, y, m nat) nat {
 	// we also multiply by x, thus adding one to the power.
 
 	w := _W - int(shift)
+	// zz and r are used to avoid allocating in mul and div as
+	// otherwise the arguments would alias.
+	var zz, r nat
 	for j := 0; j < w; j++ {
-		z = z.mul(z, z)
+		zz = zz.mul(z, z)
+		zz, z = z, zz
 
 		if v&mask != 0 {
-			z = z.mul(z, x)
+			zz = zz.mul(z, x)
+			zz, z = z, zz
 		}
 
 		if m != nil {
-			q, z = q.div(z, z, m)
+			zz, r = zz.div(r, z, m)
+			zz, r, q, z = q, z, zz, r
 		}
 
 		v <<= 1
@@ -1277,14 +1288,17 @@ func (z nat) expNN(x, y, m nat) nat {
 		v = y[i]
 
 		for j := 0; j < _W; j++ {
-			z = z.mul(z, z)
+			zz = zz.mul(z, z)
+			zz, z = z, zz
 
 			if v&mask != 0 {
-				z = z.mul(z, x)
+				zz = zz.mul(z, x)
+				zz, z = z, zz
 			}
 
 			if m != nil {
-				q, z = q.div(z, z, m)
+				zz, r = zz.div(r, z, m)
+				zz, r, q, z = q, z, zz, r
 			}
 
 			v <<= 1
@@ -1343,8 +1357,9 @@ func (n nat) probablyPrime(reps int) bool {
 	}
 
 	nm1 := nat(nil).sub(n, natOne)
-	// 1<<k * q = nm1;
-	q, k := nm1.powersOfTwoDecompose()
+	// determine q, k such that nm1 = q << k
+	k := nm1.trailingZeroBits()
+	q := nat(nil).shr(nm1, k)
 
 	nm3 := nat(nil).sub(nm1, natTwo)
 	rand := rand.New(rand.NewSource(int64(n[0])))
@@ -1360,7 +1375,7 @@ NextRandom:
 		if y.cmp(natOne) == 0 || y.cmp(nm1) == 0 {
 			continue
 		}
-		for j := 1; j < k; j++ {
+		for j := uint(1); j < k; j++ {
 			y = y.mul(y, y)
 			quotient, y = quotient.div(y, y, n)
 			if y.cmp(nm1) == 0 {
diff --git a/libgo/go/math/big/nat_test.go b/libgo/go/math/big/nat_test.go
index becde5d1719..68dd1a96d3c 100644
--- a/libgo/go/math/big/nat_test.go
+++ b/libgo/go/math/big/nat_test.go
@@ -6,6 +6,7 @@ package big
 
 import (
 	"io"
+	"runtime"
 	"strings"
 	"testing"
 )
@@ -62,6 +63,36 @@ var prodNN = []argNN{
 	{nat{0, 0, 991 * 991}, nat{0, 991}, nat{0, 991}},
 	{nat{1 * 991, 2 * 991, 3 * 991, 4 * 991}, nat{1, 2, 3, 4}, nat{991}},
 	{nat{4, 11, 20, 30, 20, 11, 4}, nat{1, 2, 3, 4}, nat{4, 3, 2, 1}},
+	// 3^100 * 3^28 = 3^128
+	{
+		natFromString("11790184577738583171520872861412518665678211592275841109096961"),
+		natFromString("515377520732011331036461129765621272702107522001"),
+		natFromString("22876792454961"),
+	},
+	// z = 111....1 (70000 digits)
+	// x = 10^(99*700) + ... + 10^1400 + 10^700 + 1
+	// y = 111....1 (700 digits, larger than Karatsuba threshold on 32-bit and 64-bit)
+	{
+		natFromString(strings.Repeat("1", 70000)),
+		natFromString("1" + strings.Repeat(strings.Repeat("0", 699)+"1", 99)),
+		natFromString(strings.Repeat("1", 700)),
+	},
+	// z = 111....1 (20000 digits)
+	// x = 10^10000 + 1
+	// y = 111....1 (10000 digits)
+	{
+		natFromString(strings.Repeat("1", 20000)),
+		natFromString("1" + strings.Repeat("0", 9999) + "1"),
+		natFromString(strings.Repeat("1", 10000)),
+	},
+}
+
+func natFromString(s string) nat {
+	x, _, err := nat(nil).scan(strings.NewReader(s), 0)
+	if err != nil {
+		panic(err)
+	}
+	return x
 }
 
 func TestSet(t *testing.T) {
@@ -135,26 +166,42 @@ func TestMulRangeN(t *testing.T) {
 	}
 }
 
-var mulArg, mulTmp nat
-
-func init() {
-	const n = 1000
-	mulArg = make(nat, n)
-	for i := 0; i < n; i++ {
-		mulArg[i] = _M
+// allocBytes returns the number of bytes allocated by invoking f. 
+func allocBytes(f func()) uint64 {
+	var stats runtime.MemStats
+	runtime.ReadMemStats(&stats)
+	t := stats.TotalAlloc
+	f()
+	runtime.ReadMemStats(&stats)
+	return stats.TotalAlloc - t
+}
+
+// TestMulUnbalanced tests that multiplying numbers of different lengths
+// does not cause deep recursion and in turn allocate too much memory.
+// Test case for issue 3807.
+func TestMulUnbalanced(t *testing.T) {
+	x := rndNat(50000)
+	y := rndNat(40)
+	allocSize := allocBytes(func() {
+		nat(nil).mul(x, y)
+	})
+	inputSize := uint64(len(x)+len(y)) * _S
+	if ratio := allocSize / uint64(inputSize); ratio > 10 {
+		t.Errorf("multiplication uses too much memory (%d > %d times the size of inputs)", allocSize, ratio)
 	}
 }
 
-func benchmarkMulLoad() {
-	for j := 1; j <= 10; j++ {
-		x := mulArg[0 : j*100]
-		mulTmp.mul(x, x)
-	}
+func rndNat(n int) nat {
+	return nat(rndV(n)).norm()
 }
 
 func BenchmarkMul(b *testing.B) {
+	mulx := rndNat(1e4)
+	muly := rndNat(1e4)
+	b.ResetTimer()
 	for i := 0; i < b.N; i++ {
-		benchmarkMulLoad()
+		var z nat
+		z.mul(mulx, muly)
 	}
 }
 
@@ -616,14 +663,23 @@ func TestModW(t *testing.T) {
 }
 
 func TestTrailingZeroBits(t *testing.T) {
-	var x Word
-	x--
-	for i := 0; i < _W; i++ {
-		if trailingZeroBits(x) != i {
-			t.Errorf("Failed at step %d: x: %x got: %d", i, x, trailingZeroBits(x))
+	x := Word(1)
+	for i := uint(0); i <= _W; i++ {
+		n := trailingZeroBits(x)
+		if n != i%_W {
+			t.Errorf("got trailingZeroBits(%#x) = %d; want %d", x, n, i%_W)
 		}
 		x <<= 1
 	}
+
+	y := nat(nil).set(natOne)
+	for i := uint(0); i <= 3*_W; i++ {
+		n := y.trailingZeroBits()
+		if n != i {
+			t.Errorf("got 0x%s.trailingZeroBits() = %d; want %d", y.string(lowercaseDigits[0:16]), n, i)
+		}
+		y = y.shl(y, 1)
+	}
 }
 
 var expNNTests = []struct {
diff --git a/libgo/go/math/big/rat.go b/libgo/go/math/big/rat.go
index 7bd83fc0fb2..eccf34e482f 100644
--- a/libgo/go/math/big/rat.go
+++ b/libgo/go/math/big/rat.go
@@ -16,8 +16,10 @@ import (
 // A Rat represents a quotient a/b of arbitrary precision.
 // The zero value for a Rat represents the value 0.
 type Rat struct {
-	a Int
-	b nat // len(b) == 0 acts like b == 1
+	// To make zero values for Rat work w/o initialization,
+	// a zero value of b (len(b) == 0) acts like b == 1.
+	// a.neg determines the sign of the Rat, b.neg is ignored.
+	a, b Int
 }
 
 // NewRat creates a new Rat with numerator a and denominator b.
@@ -36,7 +38,7 @@ func (z *Rat) SetFrac(a, b *Int) *Rat {
 		babs = nat(nil).set(babs) // make a copy
 	}
 	z.a.abs = z.a.abs.set(a.abs)
-	z.b = z.b.set(babs)
+	z.b.abs = z.b.abs.set(babs)
 	return z.norm()
 }
 
@@ -50,21 +52,21 @@ func (z *Rat) SetFrac64(a, b int64) *Rat {
 		b = -b
 		z.a.neg = !z.a.neg
 	}
-	z.b = z.b.setUint64(uint64(b))
+	z.b.abs = z.b.abs.setUint64(uint64(b))
 	return z.norm()
 }
 
 // SetInt sets z to x (by making a copy of x) and returns z.
 func (z *Rat) SetInt(x *Int) *Rat {
 	z.a.Set(x)
-	z.b = z.b.make(0)
+	z.b.abs = z.b.abs.make(0)
 	return z
 }
 
 // SetInt64 sets z to x and returns z.
 func (z *Rat) SetInt64(x int64) *Rat {
 	z.a.SetInt64(x)
-	z.b = z.b.make(0)
+	z.b.abs = z.b.abs.make(0)
 	return z
 }
 
@@ -72,7 +74,7 @@ func (z *Rat) SetInt64(x int64) *Rat {
 func (z *Rat) Set(x *Rat) *Rat {
 	if z != x {
 		z.a.Set(&x.a)
-		z.b = z.b.set(x.b)
+		z.b.Set(&x.b)
 	}
 	return z
 }
@@ -97,15 +99,15 @@ func (z *Rat) Inv(x *Rat) *Rat {
 		panic("division by zero")
 	}
 	z.Set(x)
-	a := z.b
+	a := z.b.abs
 	if len(a) == 0 {
-		a = a.setWord(1) // materialize numerator
+		a = a.set(natOne) // materialize numerator
 	}
 	b := z.a.abs
 	if b.cmp(natOne) == 0 {
 		b = b.make(0) // normalize denominator
 	}
-	z.a.abs, z.b = a, b // sign doesn't change
+	z.a.abs, z.b.abs = a, b // sign doesn't change
 	return z
 }
 
@@ -121,38 +123,26 @@ func (x *Rat) Sign() int {
 
 // IsInt returns true if the denominator of x is 1.
 func (x *Rat) IsInt() bool {
-	return len(x.b) == 0 || x.b.cmp(natOne) == 0
+	return len(x.b.abs) == 0 || x.b.abs.cmp(natOne) == 0
 }
 
 // Num returns the numerator of x; it may be <= 0.
 // The result is a reference to x's numerator; it
-// may change if a new value is assigned to x.
+// may change if a new value is assigned to x, and vice versa.
+// The sign of the numerator corresponds to the sign of x.
 func (x *Rat) Num() *Int {
 	return &x.a
 }
 
 // Denom returns the denominator of x; it is always > 0.
 // The result is a reference to x's denominator; it
-// may change if a new value is assigned to x.
+// may change if a new value is assigned to x, and vice versa.
 func (x *Rat) Denom() *Int {
-	if len(x.b) == 0 {
-		return &Int{abs: nat{1}}
+	x.b.neg = false // the result is always >= 0
+	if len(x.b.abs) == 0 {
+		x.b.abs = x.b.abs.set(natOne) // materialize denominator
 	}
-	return &Int{abs: x.b}
-}
-
-func gcd(x, y nat) nat {
-	// Euclidean algorithm.
-	var a, b nat
-	a = a.set(x)
-	b = b.set(y)
-	for len(b) != 0 {
-		var q, r nat
-		_, r = q.div(r, a, b)
-		a = b
-		b = r
-	}
-	return a
+	return &x.b
 }
 
 func (z *Rat) norm() *Rat {
@@ -160,17 +150,25 @@ func (z *Rat) norm() *Rat {
 	case len(z.a.abs) == 0:
 		// z == 0 - normalize sign and denominator
 		z.a.neg = false
-		z.b = z.b.make(0)
-	case len(z.b) == 0:
+		z.b.abs = z.b.abs.make(0)
+	case len(z.b.abs) == 0:
 		// z is normalized int - nothing to do
-	case z.b.cmp(natOne) == 0:
+	case z.b.abs.cmp(natOne) == 0:
 		// z is int - normalize denominator
-		z.b = z.b.make(0)
+		z.b.abs = z.b.abs.make(0)
 	default:
-		if f := gcd(z.a.abs, z.b); f.cmp(natOne) != 0 {
-			z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f)
-			z.b, _ = z.b.div(nil, z.b, f)
+		neg := z.a.neg
+		z.a.neg = false
+		z.b.neg = false
+		if f := NewInt(0).binaryGCD(&z.a, &z.b); f.Cmp(intOne) != 0 {
+			z.a.abs, _ = z.a.abs.div(nil, z.a.abs, f.abs)
+			z.b.abs, _ = z.b.abs.div(nil, z.b.abs, f.abs)
+			if z.b.abs.cmp(natOne) == 0 {
+				// z is int - normalize denominator
+				z.b.abs = z.b.abs.make(0)
+			}
 		}
+		z.a.neg = neg
 	}
 	return z
 }
@@ -207,31 +205,31 @@ func scaleDenom(x *Int, f nat) *Int {
 //   +1 if x >  y
 //
 func (x *Rat) Cmp(y *Rat) int {
-	return scaleDenom(&x.a, y.b).Cmp(scaleDenom(&y.a, x.b))
+	return scaleDenom(&x.a, y.b.abs).Cmp(scaleDenom(&y.a, x.b.abs))
 }
 
 // Add sets z to the sum x+y and returns z.
 func (z *Rat) Add(x, y *Rat) *Rat {
-	a1 := scaleDenom(&x.a, y.b)
-	a2 := scaleDenom(&y.a, x.b)
+	a1 := scaleDenom(&x.a, y.b.abs)
+	a2 := scaleDenom(&y.a, x.b.abs)
 	z.a.Add(a1, a2)
-	z.b = mulDenom(z.b, x.b, y.b)
+	z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
 	return z.norm()
 }
 
 // Sub sets z to the difference x-y and returns z.
 func (z *Rat) Sub(x, y *Rat) *Rat {
-	a1 := scaleDenom(&x.a, y.b)
-	a2 := scaleDenom(&y.a, x.b)
+	a1 := scaleDenom(&x.a, y.b.abs)
+	a2 := scaleDenom(&y.a, x.b.abs)
 	z.a.Sub(a1, a2)
-	z.b = mulDenom(z.b, x.b, y.b)
+	z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
 	return z.norm()
 }
 
 // Mul sets z to the product x*y and returns z.
 func (z *Rat) Mul(x, y *Rat) *Rat {
 	z.a.Mul(&x.a, &y.a)
-	z.b = mulDenom(z.b, x.b, y.b)
+	z.b.abs = mulDenom(z.b.abs, x.b.abs, y.b.abs)
 	return z.norm()
 }
 
@@ -241,10 +239,10 @@ func (z *Rat) Quo(x, y *Rat) *Rat {
 	if len(y.a.abs) == 0 {
 		panic("division by zero")
 	}
-	a := scaleDenom(&x.a, y.b)
-	b := scaleDenom(&y.a, x.b)
+	a := scaleDenom(&x.a, y.b.abs)
+	b := scaleDenom(&y.a, x.b.abs)
 	z.a.abs = a.abs
-	z.b = b.abs
+	z.b.abs = b.abs
 	z.a.neg = a.neg != b.neg
 	return z.norm()
 }
@@ -286,7 +284,7 @@ func (z *Rat) SetString(s string) (*Rat, bool) {
 		}
 		s = s[sep+1:]
 		var err error
-		if z.b, _, err = z.b.scan(strings.NewReader(s), 10); err != nil {
+		if z.b.abs, _, err = z.b.abs.scan(strings.NewReader(s), 10); err != nil {
 			return nil, false
 		}
 		return z.norm(), true
@@ -317,11 +315,11 @@ func (z *Rat) SetString(s string) (*Rat, bool) {
 	}
 	powTen := nat(nil).expNN(natTen, exp.abs, nil)
 	if exp.neg {
-		z.b = powTen
+		z.b.abs = powTen
 		z.norm()
 	} else {
 		z.a.abs = z.a.abs.mul(z.a.abs, powTen)
-		z.b = z.b.make(0)
+		z.b.abs = z.b.abs.make(0)
 	}
 
 	return z, true
@@ -330,8 +328,8 @@ func (z *Rat) SetString(s string) (*Rat, bool) {
 // String returns a string representation of z in the form "a/b" (even if b == 1).
 func (x *Rat) String() string {
 	s := "/1"
-	if len(x.b) != 0 {
-		s = "/" + x.b.decimalString()
+	if len(x.b.abs) != 0 {
+		s = "/" + x.b.abs.decimalString()
 	}
 	return x.a.String() + s
 }
@@ -355,9 +353,9 @@ func (x *Rat) FloatString(prec int) string {
 		}
 		return s
 	}
-	// x.b != 0
+	// x.b.abs != 0
 
-	q, r := nat(nil).div(nat(nil), x.a.abs, x.b)
+	q, r := nat(nil).div(nat(nil), x.a.abs, x.b.abs)
 
 	p := natOne
 	if prec > 0 {
@@ -365,11 +363,11 @@ func (x *Rat) FloatString(prec int) string {
 	}
 
 	r = r.mul(r, p)
-	r, r2 := r.div(nat(nil), r, x.b)
+	r, r2 := r.div(nat(nil), r, x.b.abs)
 
 	// see if we need to round up
 	r2 = r2.add(r2, r2)
-	if x.b.cmp(r2) <= 0 {
+	if x.b.abs.cmp(r2) <= 0 {
 		r = r.add(r, natOne)
 		if r.cmp(p) >= 0 {
 			q = nat(nil).add(q, natOne)
@@ -396,8 +394,8 @@ const ratGobVersion byte = 1
 
 // GobEncode implements the gob.GobEncoder interface.
 func (x *Rat) GobEncode() ([]byte, error) {
-	buf := make([]byte, 1+4+(len(x.a.abs)+len(x.b))*_S) // extra bytes for version and sign bit (1), and numerator length (4)
-	i := x.b.bytes(buf)
+	buf := make([]byte, 1+4+(len(x.a.abs)+len(x.b.abs))*_S) // extra bytes for version and sign bit (1), and numerator length (4)
+	i := x.b.abs.bytes(buf)
 	j := x.a.abs.bytes(buf[0:i])
 	n := i - j
 	if int(uint32(n)) != n {
@@ -427,6 +425,6 @@ func (z *Rat) GobDecode(buf []byte) error {
 	i := j + binary.BigEndian.Uint32(buf[j-4:j])
 	z.a.neg = b&1 != 0
 	z.a.abs = z.a.abs.setBytes(buf[j:i])
-	z.b = z.b.setBytes(buf[i:])
+	z.b.abs = z.b.abs.setBytes(buf[i:])
 	return nil
 }
diff --git a/libgo/go/math/big/rat_test.go b/libgo/go/math/big/rat_test.go
index f7f31ae1a20..7c634233ff8 100644
--- a/libgo/go/math/big/rat_test.go
+++ b/libgo/go/math/big/rat_test.go
@@ -387,30 +387,19 @@ func TestRatGobEncoding(t *testing.T) {
 	var medium bytes.Buffer
 	enc := gob.NewEncoder(&medium)
 	dec := gob.NewDecoder(&medium)
-	for i, test := range gobEncodingTests {
-		for j := 0; j < 4; j++ {
-			medium.Reset() // empty buffer for each test case (in case of failures)
-			stest := test
-			if j&1 != 0 {
-				// negative numbers
-				stest = "-" + test
-			}
-			if j%2 != 0 {
-				// fractions
-				stest = stest + "." + test
-			}
-			var tx Rat
-			tx.SetString(stest)
-			if err := enc.Encode(&tx); err != nil {
-				t.Errorf("#%d%c: encoding failed: %s", i, 'a'+j, err)
-			}
-			var rx Rat
-			if err := dec.Decode(&rx); err != nil {
-				t.Errorf("#%d%c: decoding failed: %s", i, 'a'+j, err)
-			}
-			if rx.Cmp(&tx) != 0 {
-				t.Errorf("#%d%c: transmission failed: got %s want %s", i, 'a'+j, &rx, &tx)
-			}
+	for _, test := range encodingTests {
+		medium.Reset() // empty buffer for each test case (in case of failures)
+		var tx Rat
+		tx.SetString(test + ".14159265")
+		if err := enc.Encode(&tx); err != nil {
+			t.Errorf("encoding of %s failed: %s", &tx, err)
+		}
+		var rx Rat
+		if err := dec.Decode(&rx); err != nil {
+			t.Errorf("decoding of %s failed: %s", &tx, err)
+		}
+		if rx.Cmp(&tx) != 0 {
+			t.Errorf("transmission of %s failed: got %s want %s", &tx, &rx, &tx)
 		}
 	}
 }
@@ -454,3 +443,56 @@ func TestIssue2379(t *testing.T) {
 		t.Errorf("5) got %s want %s", x, q)
 	}
 }
+
+func TestIssue3521(t *testing.T) {
+	a := new(Int)
+	b := new(Int)
+	a.SetString("64375784358435883458348587", 0)
+	b.SetString("4789759874531", 0)
+
+	// 0) a raw zero value has 1 as denominator
+	zero := new(Rat)
+	one := NewInt(1)
+	if zero.Denom().Cmp(one) != 0 {
+		t.Errorf("0) got %s want %s", zero.Denom(), one)
+	}
+
+	// 1a) a zero value remains zero independent of denominator
+	x := new(Rat)
+	x.Denom().Set(new(Int).Neg(b))
+	if x.Cmp(zero) != 0 {
+		t.Errorf("1a) got %s want %s", x, zero)
+	}
+
+	// 1b) a zero value may have a denominator != 0 and != 1
+	x.Num().Set(a)
+	qab := new(Rat).SetFrac(a, b)
+	if x.Cmp(qab) != 0 {
+		t.Errorf("1b) got %s want %s", x, qab)
+	}
+
+	// 2a) an integral value becomes a fraction depending on denominator
+	x.SetFrac64(10, 2)
+	x.Denom().SetInt64(3)
+	q53 := NewRat(5, 3)
+	if x.Cmp(q53) != 0 {
+		t.Errorf("2a) got %s want %s", x, q53)
+	}
+
+	// 2b) an integral value becomes a fraction depending on denominator
+	x = NewRat(10, 2)
+	x.Denom().SetInt64(3)
+	if x.Cmp(q53) != 0 {
+		t.Errorf("2b) got %s want %s", x, q53)
+	}
+
+	// 3) changing the numerator/denominator of a Rat changes the Rat
+	x.SetFrac(a, b)
+	a = x.Num()
+	b = x.Denom()
+	a.SetInt64(5)
+	b.SetInt64(3)
+	if x.Cmp(q53) != 0 {
+		t.Errorf("3) got %s want %s", x, q53)
+	}
+}
diff --git a/libgo/go/math/cbrt.go b/libgo/go/math/cbrt.go
index 8c43f0afbc8..272e3092310 100644
--- a/libgo/go/math/cbrt.go
+++ b/libgo/go/math/cbrt.go
@@ -12,7 +12,7 @@ package math
 	(http://www.jstor.org/stable/2006387?seq=9, accessed 11-Feb-2010)
 */
 
-// Cbrt returns the cube root of its argument.
+// Cbrt returns the cube root of x.
 //
 // Special cases are:
 //	Cbrt(±0) = ±0
diff --git a/libgo/go/math/copysign.go b/libgo/go/math/copysign.go
index ee65456a1ce..719c64b9eba 100644
--- a/libgo/go/math/copysign.go
+++ b/libgo/go/math/copysign.go
@@ -4,7 +4,7 @@
 
 package math
 
-// Copysign(x, y) returns a value with the magnitude
+// Copysign returns a value with the magnitude
 // of x and the sign of y.
 func Copysign(x, y float64) float64 {
 	const sign = 1 << 63
diff --git a/libgo/go/math/erf.go b/libgo/go/math/erf.go
index c6f32bdbe28..4cd80f80c3b 100644
--- a/libgo/go/math/erf.go
+++ b/libgo/go/math/erf.go
@@ -179,7 +179,7 @@ const (
 	sb7 = -2.24409524465858183362e+01 // 0xC03670E242712D62
 )
 
-// Erf(x) returns the error function of x.
+// Erf returns the error function of x.
 //
 // Special cases are:
 //	Erf(+Inf) = 1
@@ -256,7 +256,7 @@ func Erf(x float64) float64 {
 	return 1 - r/x
 }
 
-// Erfc(x) returns the complementary error function of x.
+// Erfc returns the complementary error function of x.
 //
 // Special cases are:
 //	Erfc(+Inf) = 0
diff --git a/libgo/go/math/gamma.go b/libgo/go/math/gamma.go
index 7c6f421bad1..164f54f332b 100644
--- a/libgo/go/math/gamma.go
+++ b/libgo/go/math/gamma.go
@@ -110,19 +110,26 @@ func stirling(x float64) float64 {
 	return y
 }
 
-// Gamma(x) returns the Gamma function of x.
+// Gamma returns the Gamma function of x.
 //
 // Special cases are:
-//	Gamma(±Inf) = ±Inf
+//	Gamma(+Inf) = +Inf
+//	Gamma(+0) = +Inf
+//	Gamma(-0) = -Inf
+//	Gamma(x) = NaN for integer x < 0
+//	Gamma(-Inf) = NaN
 //	Gamma(NaN) = NaN
-// Large values overflow to +Inf.
-// Zero and negative integer arguments return ±Inf.
 func Gamma(x float64) float64 {
 	const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620
 	// special cases
 	switch {
-	case IsInf(x, -1) || IsNaN(x):
-		return x
+	case isNegInt(x) || IsInf(x, -1) || IsNaN(x):
+		return NaN()
+	case x == 0:
+		if Signbit(x) {
+			return Inf(-1)
+		}
+		return Inf(1)
 	case x < -170.5674972726612 || x > 171.61447887182298:
 		return Inf(1)
 	}
@@ -185,3 +192,11 @@ small:
 	}
 	return z / ((1 + Euler*x) * x)
 }
+
+func isNegInt(x float64) bool {
+	if x < 0 {
+		_, xf := Modf(x)
+		return xf == 0
+	}
+	return false
+}
diff --git a/libgo/go/math/hypot.go b/libgo/go/math/hypot.go
index 57d8e343720..a6fa84c7655 100644
--- a/libgo/go/math/hypot.go
+++ b/libgo/go/math/hypot.go
@@ -8,7 +8,7 @@ package math
 	Hypot -- sqrt(p*p + q*q), but overflows only if the result does.
 */
 
-// Hypot computes Sqrt(p*p + q*q), taking care to avoid
+// Hypot returns Sqrt(p*p + q*q), taking care to avoid
 // unnecessary overflow and underflow.
 //
 // Special cases are:
diff --git a/libgo/go/math/logb.go b/libgo/go/math/logb.go
index d32f9f1000c..f2769d4fd75 100644
--- a/libgo/go/math/logb.go
+++ b/libgo/go/math/logb.go
@@ -4,7 +4,7 @@
 
 package math
 
-// Logb(x) returns the binary exponent of x.
+// Logb returns the binary exponent of x.
 //
 // Special cases are:
 //	Logb(±Inf) = +Inf
@@ -23,7 +23,7 @@ func Logb(x float64) float64 {
 	return float64(ilogb(x))
 }
 
-// Ilogb(x) returns the binary exponent of x as an integer.
+// Ilogb returns the binary exponent of x as an integer.
 //
 // Special cases are:
 //	Ilogb(±Inf) = MaxInt32
diff --git a/libgo/go/math/rand/rand_test.go b/libgo/go/math/rand/rand_test.go
index bbd44e3f8b1..4d3abdb606c 100644
--- a/libgo/go/math/rand/rand_test.go
+++ b/libgo/go/math/rand/rand_test.go
@@ -57,16 +57,13 @@ func (this *statsResults) checkSimilarDistribution(expected *statsResults) error
 
 func getStatsResults(samples []float64) *statsResults {
 	res := new(statsResults)
-	var sum float64
-	for i := range samples {
-		sum += samples[i]
+	var sum, squaresum float64
+	for _, s := range samples {
+		sum += s
+		squaresum += s * s
 	}
 	res.mean = sum / float64(len(samples))
-	var devsum float64
-	for i := range samples {
-		devsum += math.Pow(samples[i]-res.mean, 2)
-	}
-	res.stddev = math.Sqrt(devsum / float64(len(samples)))
+	res.stddev = math.Sqrt(squaresum/float64(len(samples)) - res.mean*res.mean)
 	return res
 }
 
diff --git a/libgo/go/math/sincos.go b/libgo/go/math/sincos.go
index 75e6e7541e6..b3a2f8a8619 100644
--- a/libgo/go/math/sincos.go
+++ b/libgo/go/math/sincos.go
@@ -6,7 +6,7 @@ package math
 
 // Coefficients _sin[] and _cos[] are found in pkg/math/sin.go.
 
-// Sincos(x) returns Sin(x), Cos(x).
+// Sincos returns Sin(x), Cos(x).
 //
 // Special cases are:
 //	Sincos(±0) = ±0, 1
diff --git a/libgo/go/math/tanh.go b/libgo/go/math/tanh.go
index 03a641b4da0..7305be66c7f 100644
--- a/libgo/go/math/tanh.go
+++ b/libgo/go/math/tanh.go
@@ -4,12 +4,66 @@
 
 package math
 
-/*
-	Floating-point hyperbolic tangent.
+// The original C code, the long comment, and the constants
+// below were from http://netlib.sandia.gov/cephes/cmath/sin.c,
+// available from http://www.netlib.org/cephes/cmath.tgz.
+// The go code is a simplified version of the original C.
+//      tanh.c
+//
+//      Hyperbolic tangent
+//
+// SYNOPSIS:
+//
+// double x, y, tanh();
+//
+// y = tanh( x );
+//
+// DESCRIPTION:
+//
+// Returns hyperbolic tangent of argument in the range MINLOG to MAXLOG.
+//      MAXLOG = 8.8029691931113054295988e+01 = log(2**127)
+//      MINLOG = -8.872283911167299960540e+01 = log(2**-128)
+//
+// A rational function is used for |x| < 0.625.  The form
+// x + x**3 P(x)/Q(x) of Cody & Waite is employed.
+// Otherwise,
+//      tanh(x) = sinh(x)/cosh(x) = 1  -  2/(exp(2x) + 1).
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain     # trials      peak         rms
+//    IEEE      -2,2        30000       2.5e-16     5.8e-17
+//
+// Cephes Math Library Release 2.8:  June, 2000
+// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+//
+// The readme file at http://netlib.sandia.gov/cephes/ says:
+//    Some software in this archive may be from the book _Methods and
+// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
+// International, 1989) or from the Cephes Mathematical Library, a
+// commercial product. In either event, it is copyrighted by the author.
+// What you see here may be used freely but it comes with no support or
+// guarantee.
+//
+//   The two known misprints in the book are repaired here in the
+// source listings for the gamma function and the incomplete beta
+// integral.
+//
+//   Stephen L. Moshier
+//   moshier@na-net.ornl.gov
+//
 
-	Sinh and Cosh are called except for large arguments, which
-	would cause overflow improperly.
-*/
+var tanhP = [...]float64{
+	-9.64399179425052238628E-1,
+	-9.92877231001918586564E1,
+	-1.61468768441708447952E3,
+}
+var tanhQ = [...]float64{
+	1.12811678491632931402E2,
+	2.23548839060100448583E3,
+	4.84406305325125486048E3,
+}
 
 // Tanh computes the hyperbolic tangent of x.
 //
@@ -18,15 +72,26 @@ package math
 //	Tanh(±Inf) = ±1
 //	Tanh(NaN) = NaN
 func Tanh(x float64) float64 {
-	if x < 0 {
-		x = -x
-		if x > 21 {
+	const MAXLOG = 8.8029691931113054295988e+01 // log(2**127)
+	z := Abs(x)
+	switch {
+	case z > 0.5*MAXLOG:
+		if x < 0 {
 			return -1
 		}
-		return -Sinh(x) / Cosh(x)
-	}
-	if x > 21 {
 		return 1
+	case z >= 0.625:
+		s := Exp(2 * z)
+		z = 1 - 2/(s+1)
+		if x < 0 {
+			z = -z
+		}
+	default:
+		if x == 0 {
+			return x
+		}
+		s := x * x
+		z = x + x*s*((tanhP[0]*s+tanhP[1])*s+tanhP[2])/(((s+tanhQ[0])*s+tanhQ[1])*s+tanhQ[2])
 	}
-	return Sinh(x) / Cosh(x)
+	return z
 }