math: special cases for Pow

R=rsc
CC=golang-dev
https://golang.org/cl/176064

2 files changed, 202 insertions, 5 deletions

diff --git a/src/pkg/math/all_test.go b/src/pkg/math/all_test.go
index 6feddac932..dc6177dad6 100644
--- a/src/pkg/math/all_test.go
+++ b/src/pkg/math/all_test.go
@@ -154,6 +154,126 @@ var tanh = []float64{
 	9.4936501296239700e-01,
 	-9.9999994291374019e-01,
 }
+var vfsin = []float64{
+	NaN(),
+	Inf(-1),
+	0,
+	Inf(1),
+}
+var vfasin = []float64{
+	NaN(),
+	-Pi,
+	0,
+	Pi,
+}
+var vf1 = []float64{
+	NaN(),
+	Inf(-1),
+	-Pi,
+	-1,
+	0,
+	1,
+	Pi,
+	Inf(1),
+}
+var vfhypot = [][2]float64{
+	[2]float64{Inf(-1), 1},
+	[2]float64{Inf(1), 1},
+	[2]float64{1, Inf(-1)},
+	[2]float64{1, Inf(1)},
+	[2]float64{NaN(), Inf(-1)},
+	[2]float64{NaN(), Inf(1)},
+	[2]float64{1, NaN()},
+	[2]float64{NaN(), 1},
+}
+var vf2 = [][2]float64{
+	[2]float64{-Pi, Pi},
+	[2]float64{-Pi, -Pi},
+	[2]float64{Inf(-1), 3},
+	[2]float64{Inf(-1), Pi},
+	[2]float64{Inf(-1), -3},
+	[2]float64{Inf(-1), -Pi},
+	[2]float64{Inf(1), Pi},
+	[2]float64{0, -Pi},
+	[2]float64{Inf(1), -Pi},
+	[2]float64{0, Pi},
+	[2]float64{-1, Inf(-1)},
+	[2]float64{-1, Inf(1)},
+	[2]float64{1, Inf(-1)},
+	[2]float64{1, Inf(1)},
+	[2]float64{-1 / 2, Inf(1)},
+	[2]float64{1 / 2, Inf(1)},
+	[2]float64{-Pi, Inf(-1)},
+	[2]float64{Pi, Inf(-1)},
+	[2]float64{-1 / 2, Inf(-1)},
+	[2]float64{1 / 2, Inf(-1)},
+	[2]float64{-Pi, Inf(1)},
+	[2]float64{Pi, Inf(1)},
+	[2]float64{NaN(), -Pi},
+	[2]float64{NaN(), Pi},
+	[2]float64{Inf(-1), NaN()},
+	[2]float64{-Pi, NaN()},
+	[2]float64{0, NaN()},
+	[2]float64{Pi, NaN()},
+	[2]float64{Inf(1), NaN()},
+	[2]float64{NaN(), NaN()},
+	[2]float64{Inf(-1), 1},
+	[2]float64{-Pi, 1},
+	[2]float64{0, 1},
+	[2]float64{Pi, 1},
+	[2]float64{Inf(1), 1},
+	[2]float64{NaN(), 1},
+	[2]float64{Inf(-1), 0},
+	[2]float64{-Pi, 0},
+	[2]float64{0, 0},
+	[2]float64{Pi, 0},
+	[2]float64{Inf(1), 0},
+	[2]float64{NaN(), 0},
+}
+var pow2 = []float64{
+	NaN(),
+	NaN(),
+	Inf(-1),
+	Inf(1),
+	0,
+	0,
+	Inf(1),
+	Inf(1),
+	0,
+	0,
+	NaN(),
+	NaN(),
+	NaN(),
+	NaN(),
+	0,
+	0,
+	0,
+	0,
+	Inf(1),
+	Inf(1),
+	Inf(1),
+	Inf(1),
+	NaN(),
+	NaN(),
+	NaN(),
+	NaN(),
+	NaN(),
+	NaN(),
+	NaN(),
+	NaN(),
+	Inf(-1),
+	-Pi,
+	0,
+	Pi,
+	Inf(1),
+	NaN(),
+	1,
+	1,
+	1,
+	1,
+	1,
+	1,
+}
 
 func tolerance(a, b, e float64) bool {
 	d := a - b
@@ -172,6 +292,19 @@ func tolerance(a, b, e float64) bool {
 func kindaclose(a, b float64) bool { return tolerance(a, b, 1e-8) }
 func close(a, b float64) bool      { return tolerance(a, b, 1e-14) }
 func veryclose(a, b float64) bool  { return tolerance(a, b, 4e-16) }
+func alike(a, b float64) bool {
+	switch {
+	case IsNaN(a) && IsNaN(b):
+		return true
+	case IsInf(a, 1) && IsInf(b, 1):
+		return true
+	case IsInf(a, -1) && IsInf(b, -1):
+		return true
+	case a == b:
+		return true
+	}
+	return false
+}
 
 func TestAsin(t *testing.T) {
 	for i := 0; i < len(vf); i++ {
@@ -223,6 +356,11 @@ func TestPow(t *testing.T) {
 			t.Errorf("Pow(10, %.17g) = %.17g, want %.17g\n", vf[i], f, pow[i])
 		}
 	}
+	for i := 0; i < len(vf2); i++ {
+		if f := Pow(vf2[i][0], vf2[i][1]); !alike(pow2[i], f) {
+			t.Errorf("Pow(%.17g, %.17g) = %.17g, want %.17g\n", vf2[i][0], vf2[i][1], f, pow2[i])
+		}
+	}
 }
 
 func TestSin(t *testing.T) {
@@ -336,3 +474,17 @@ func TestFloatMinMax(t *testing.T) {
 		}
 	}
 }
+
+// Benchmarks
+
+func BenchmarkPowInt(b *testing.B) {
+	for i := 0; i < b.N; i++ {
+		Pow(2, 2)
+	}
+}
+
+func BenchmarkPowFrac(b *testing.B) {
+	for i := 0; i < b.N; i++ {
+		Pow(2.5, 1.5)
+	}
+}
diff --git a/src/pkg/math/pow.go b/src/pkg/math/pow.go
index 42a968d395..ab8bdb60c1 100644
--- a/src/pkg/math/pow.go
+++ b/src/pkg/math/pow.go
@@ -4,23 +4,68 @@
 
 package math
 
+func isOddInt(x float64) bool {
+	xi, xf := Modf(x)
+	return xf == 0 && int64(xi)&1 == 1
+}
 
 // Pow returns x**y, the base-x exponential of y.
 func Pow(x, y float64) float64 {
-	// TODO: x or y NaN, ±Inf, maybe ±0.
+	// TODO:  maybe ±0.
+	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
+	// when compiler does it for us
 	switch {
 	case y == 0:
 		return 1
 	case y == 1:
 		return x
-	case x == 0 && y > 0:
-		return 0
-	case x == 0 && y < 0:
-		return Inf(1)
 	case y == 0.5:
 		return Sqrt(x)
 	case y == -0.5:
 		return 1 / Sqrt(x)
+	case x != x || y != y: // IsNaN(x) || IsNaN(y):
+		return NaN()
+	case x == 0:
+		switch {
+		case y < 0:
+			return Inf(1)
+		case y > 0:
+			return 0
+		}
+	case y > MaxFloat64 || y < -MaxFloat64: // IsInf(y, 0):
+		switch {
+		case Fabs(x) == 1:
+			return NaN()
+		case Fabs(x) < 1:
+			switch {
+			case IsInf(y, -1):
+				return Inf(1)
+			case IsInf(y, 1):
+				return 0
+			}
+		case Fabs(x) > 1:
+			switch {
+			case IsInf(y, -1):
+				return 0
+			case IsInf(y, 1):
+				return Inf(1)
+			}
+		}
+	case x > MaxFloat64 || x < -MaxFloat64: // IsInf(x, 0):
+		switch {
+		case y < 0:
+			return 0
+		case y > 0:
+			switch {
+			case IsInf(x, -1):
+				if isOddInt(y) {
+					return Inf(-1)
+				}
+				return Inf(1)
+			case IsInf(x, 1):
+				return Inf(1)
+			}
+		}
 	}
 
 	absy := y