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
http://codereview.appspot.com/229049

Committer: Russ Cox <rsc@golang.org>

5 files changed, 81 insertions, 24 deletions

diff --git a/src/pkg/math/Makefile b/src/pkg/math/Makefile
index 6650482a7..3b82a786b 100644
--- 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 627949971..8cb575659 100644
--- 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 31924165e..ecd115d9e 100644
--- 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
index 000000000..27f335ba2
--- /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
index 000000000..85ce1d404
--- /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) }