math: faster hypot

Use hardware sqrt for faster hypot; preserve software-only
hypot as hypotGo (like sqrtGo); enable benchmarking of
hypotGo.

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

diff --git a/src/pkg/math/Makefile b/src/pkg/math/Makefile
index 6650482a7e708567202bdc91f8dce4e7ca5daca7..3b82a786b3e7cce7934c61c9d917fa382ff005d8 100644 (file)
--- a/src/pkg/math/Makefile
+++ b/src/pkg/math/Makefile
@@ -54,6 +54,7 @@ ALLGOFILES=\
        fmod.go\
        frexp.go\
        hypot.go\
+       hypot_port.go\
        logb.go\
        lgamma.go\
        ldexp.go\

diff --git a/src/pkg/math/all_test.go b/src/pkg/math/all_test.go
index 6279499713a1b79ed79b08831268c92db201fa9f..8cb575659af304c1c13a1fba8d22a7bb53552617 100644 (file)
--- a/src/pkg/math/all_test.go
+++ b/src/pkg/math/all_test.go
@@ -1754,6 +1754,12 @@ func BenchmarkHypot(b *testing.B) {
        }
 }
 
+func BenchmarkHypotGo(b *testing.B) {
+       for i := 0; i < b.N; i++ {
+               HypotGo(3, 4)
+       }
+}
+
 func BenchmarkIlogb(b *testing.B) {
        for i := 0; i < b.N; i++ {
                Ilogb(.5)

diff --git a/src/pkg/math/hypot.go b/src/pkg/math/hypot.go
index 31924165e7e71f2ae127669fbfd1518a825995a4..ecd115d9ef8bc4762d3aace8d6f281423dd5a7ac 100644 (file)
--- a/src/pkg/math/hypot.go
+++ b/src/pkg/math/hypot.go
@@ -1,4 +1,4 @@
-// Copyright 2009-2010 The Go Authors. All rights reserved.
+// Copyright 2010 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.
 
@@ -6,11 +6,6 @@ package math
 
 /*
        Hypot -- sqrt(p*p + q*q), but overflows only if the result does.
-       See:
-               Cleve Moler and Donald Morrison,
-               Replacing Square Roots by Pythagorean Sums
-               IBM Journal of Research and Development,
-               Vol. 27, Number 6, pp. 577-581, Nov. 1983
 */
 
 // Hypot computes Sqrt(p*p + q*q), taking care to avoid
@@ -35,29 +30,12 @@ func Hypot(p, q float64) float64 {
        if q < 0 {
                q = -q
        }
-
        if p < q {
                p, q = q, p
        }
-
        if p == 0 {
                return 0
        }
-
-       pfac := p
        q = q / p
-       r := q
-       p = 1
-       for {
-               r = r * r
-               s := r + 4
-               if s == 4 {
-                       return p * pfac
-               }
-               r = r / s
-               p = p + 2*r*p
-               q = q * r
-               r = q / p
-       }
-       panic("unreachable")
+       return p * Sqrt(1+q*q)
 }

diff --git a/src/pkg/math/hypot_port.go b/src/pkg/math/hypot_port.go
new file mode 100644 (file)
index 0000000..27f335b
--- /dev/null
+++ b/src/pkg/math/hypot_port.go
@@ -0,0 +1,63 @@
+// Copyright 2009-2010 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 math
+
+/*
+       Hypot -- sqrt(p*p + q*q), but overflows only if the result does.
+       See:
+               Cleve Moler and Donald Morrison,
+               Replacing Square Roots by Pythagorean Sums
+               IBM Journal of Research and Development,
+               Vol. 27, Number 6, pp. 577-581, Nov. 1983
+*/
+
+// Hypot computes Sqrt(p*p + q*q), taking care to avoid
+// unnecessary overflow and underflow.
+//
+// Special cases are:
+//     Hypot(p, q) = +Inf if p or q is infinite
+//     Hypot(p, q) = NaN if p or q is NaN
+func hypotGo(p, q float64) float64 {
+       // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+       // when compiler does it for us
+       // special cases
+       switch {
+       case p < -MaxFloat64 || p > MaxFloat64 || q < -MaxFloat64 || q > MaxFloat64: // IsInf(p, 0) || IsInf(q, 0):
+               return Inf(1)
+       case p != p || q != q: // IsNaN(p) || IsNaN(q):
+               return NaN()
+       }
+       if p < 0 {
+               p = -p
+       }
+       if q < 0 {
+               q = -q
+       }
+
+       if p < q {
+               p, q = q, p
+       }
+
+       if p == 0 {
+               return 0
+       }
+
+       pfac := p
+       q = q / p
+       r := q
+       p = 1
+       for {
+               r = r * r
+               s := r + 4
+               if s == 4 {
+                       return p * pfac
+               }
+               r = r / s
+               p = p + 2*r*p
+               q = q * r
+               r = q / p
+       }
+       panic("unreachable")
+}

diff --git a/src/pkg/math/hypot_test.go b/src/pkg/math/hypot_test.go
new file mode 100644 (file)
index 0000000..85ce1d4
--- /dev/null
+++ b/src/pkg/math/hypot_test.go
@@ -0,0 +1,9 @@
+// Copyright 2010 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 math
+
+// Make hypotGo available for testing.
+
+func HypotGo(x, y float64) float64 { return hypotGo(x, y) }