math: special cases for Pow

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

diff --git a/src/pkg/math/all_test.go b/src/pkg/math/all_test.go
index 6feddac932cce6cd4486814f4c35c8fe16dc7e19..dc6177dad61a9d24e580345ec8d8d8921ffc89e6 100644 (file)
--- 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 42a968d395b38c492933d9fa3fa6c21d02f1a85b..ab8bdb60c19cdb5bc4c6311e3babbe4a7bea514f 100644 (file)
--- 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