math: Sqrt using 386 FPU.

Note: sqrt_decl.go already in src/pkg/math/.

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

5 files changed, 162 insertions, 129 deletions

diff --git a/src/pkg/math/Makefile b/src/pkg/math/Makefile
index f30f38fafe..7a3808976a 100644
--- a/src/pkg/math/Makefile
+++ b/src/pkg/math/Makefile
@@ -9,6 +9,9 @@ TARG=math
 OFILES_amd64=\
 	sqrt_amd64.$O\
 
+OFILES_386=\
+	sqrt_386.$O\
+
 OFILES=\
 	$(OFILES_$(GOARCH))
 
@@ -29,6 +32,7 @@ ALLGOFILES=\
 	sin.go\
 	sinh.go\
 	sqrt.go\
+	sqrt_port.go\
 	tan.go\
 	tanh.go\
 	unsafe.go\
diff --git a/src/pkg/math/all_test.go b/src/pkg/math/all_test.go
index 58728801b4..04c273322b 100644
--- a/src/pkg/math/all_test.go
+++ b/src/pkg/math/all_test.go
@@ -529,3 +529,9 @@ func BenchmarkAcos(b *testing.B) {
 		Acos(.5)
 	}
 }
+
+func BenchmarkSqrt(b *testing.B) {
+	for i := 0; i < b.N; i++ {
+		Sqrt(10)
+	}
+}
diff --git a/src/pkg/math/sqrt.go b/src/pkg/math/sqrt.go
index a3a3119fed..f12e48734f 100644
--- a/src/pkg/math/sqrt.go
+++ b/src/pkg/math/sqrt.go
@@ -4,84 +4,6 @@
 
 package math
 
-// The original C code and the long comment below are
-// from FreeBSD's /usr/src/lib/msun/src/e_sqrt.c and
-// came with this notice.  The go code is a simplified
-// version of the original C.
-//
-// ====================================================
-// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
-//
-// Developed at SunPro, a Sun Microsystems, Inc. business.
-// Permission to use, copy, modify, and distribute this
-// software is freely granted, provided that this notice
-// is preserved.
-// ====================================================
-//
-// __ieee754_sqrt(x)
-// Return correctly rounded sqrt.
-//           -----------------------------------------
-//           | Use the hardware sqrt if you have one |
-//           -----------------------------------------
-// Method:
-//   Bit by bit method using integer arithmetic. (Slow, but portable)
-//   1. Normalization
-//      Scale x to y in [1,4) with even powers of 2:
-//      find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
-//              sqrt(x) = 2^k * sqrt(y)
-//   2. Bit by bit computation
-//      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
-//           i                                                   0
-//                                     i+1         2
-//          s  = 2*q , and      y  =  2   * ( y - q  ).          (1)
-//           i      i            i                 i
-//
-//      To compute q    from q , one checks whether
-//                  i+1       i
-//
-//                            -(i+1) 2
-//                      (q + 2      )  <= y.                     (2)
-//                        i
-//                                                            -(i+1)
-//      If (2) is false, then q   = q ; otherwise q   = q  + 2      .
-//                             i+1   i             i+1   i
-//
-//      With some algebric manipulation, it is not difficult to see
-//      that (2) is equivalent to
-//                             -(i+1)
-//                      s  +  2       <= y                       (3)
-//                       i                i
-//
-//      The advantage of (3) is that s  and y  can be computed by
-//                                    i      i
-//      the following recurrence formula:
-//          if (3) is false
-//
-//          s     =  s  ,       y    = y   ;                     (4)
-//           i+1      i          i+1    i
-//
-//      otherwise,
-//                         -i                      -(i+1)
-//          s     =  s  + 2  ,  y    = y  -  s  - 2              (5)
-//           i+1      i          i+1    i     i
-//
-//      One may easily use induction to prove (4) and (5).
-//      Note. Since the left hand side of (3) contain only i+2 bits,
-//            it does not necessary to do a full (53-bit) comparison
-//            in (3).
-//   3. Final rounding
-//      After generating the 53 bits result, we compute one more bit.
-//      Together with the remainder, we can decide whether the
-//      result is exact, bigger than 1/2ulp, or less than 1/2ulp
-//      (it will never equal to 1/2ulp).
-//      The rounding mode can be detected by checking whether
-//      huge + tiny is equal to huge, and whether huge - tiny is
-//      equal to huge for some floating point number "huge" and "tiny".
-//
-//
-// Notes:  Rounding mode detection omitted.  The constants "mask", "shift",
-// and "bias" are found in src/pkg/math/bits.go
-
 // Sqrt returns the square root of x.
 //
 // Special cases are:
@@ -89,54 +11,4 @@ package math
 //	Sqrt(0) = 0
 //	Sqrt(x < 0) = NaN
 //	Sqrt(NaN) = NaN
-func Sqrt(x float64) float64 {
-	// special cases
-	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
-	// when compiler does it for us
-	switch {
-	case x != x || x > MaxFloat64: // IsNaN(x) || IsInf(x, 1):
-		return x
-	case x == 0:
-		return 0
-	case x < 0:
-		return NaN()
-	}
-	ix := Float64bits(x)
-	// normalize x
-	exp := int((ix >> shift) & mask)
-	if exp == 0 { // subnormal x
-		for ix&1<<shift == 0 {
-			ix <<= 1
-			exp--
-		}
-		exp++
-	}
-	exp -= bias + 1 // unbias exponent
-	ix &^= mask << shift
-	ix |= 1 << shift
-	if exp&1 == 1 { // odd exp, double x to make it even
-		ix <<= 1
-	}
-	exp >>= 1 // exp = exp/2, exponent of square root
-	// generate sqrt(x) bit by bit
-	ix <<= 1
-	var q, s uint64               // q = sqrt(x)
-	r := uint64(1 << (shift + 1)) // r = moving bit from MSB to LSB
-	for r != 0 {
-		t := s + r
-		if t <= ix {
-			s = t + r
-			ix -= t
-			q += r
-		}
-		ix <<= 1
-		r >>= 1
-	}
-	// final rounding
-	if ix != 0 { // remainder, result not exact
-		q += q & 1 // round according to extra bit
-	}
-	ix = q>>1 + 0x3fe0000000000000 // q/2 + 0.5
-	ix += uint64(exp) << shift
-	return Float64frombits(ix)
-}
+func Sqrt(x float64) float64 { return sqrtGo(x) }
diff --git a/src/pkg/math/sqrt_386.s b/src/pkg/math/sqrt_386.s
new file mode 100644
index 0000000000..c3bad21280
--- /dev/null
+++ b/src/pkg/math/sqrt_386.s
@@ -0,0 +1,10 @@
+// Copyright 2009 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.
+
+// func Sqrt(x float64) float64	
+TEXT math·Sqrt(SB),7,$0
+	FMOVD   x+0(FP),F0
+	FSQRT
+	FMOVDP  F0,r+8(FP)
+	RET
diff --git a/src/pkg/math/sqrt_port.go b/src/pkg/math/sqrt_port.go
new file mode 100644
index 0000000000..feccbc6199
--- /dev/null
+++ b/src/pkg/math/sqrt_port.go
@@ -0,0 +1,141 @@
+// Copyright 2009 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
+
+// The original C code and the long comment below are
+// from FreeBSD's /usr/src/lib/msun/src/e_sqrt.c and
+// came with this notice.  The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+// __ieee754_sqrt(x)
+// Return correctly rounded sqrt.
+//           -----------------------------------------
+//           | Use the hardware sqrt if you have one |
+//           -----------------------------------------
+// Method:
+//   Bit by bit method using integer arithmetic. (Slow, but portable)
+//   1. Normalization
+//      Scale x to y in [1,4) with even powers of 2:
+//      find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
+//              sqrt(x) = 2^k * sqrt(y)
+//   2. Bit by bit computation
+//      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
+//           i                                                   0
+//                                     i+1         2
+//          s  = 2*q , and      y  =  2   * ( y - q  ).          (1)
+//           i      i            i                 i
+//
+//      To compute q    from q , one checks whether
+//                  i+1       i
+//
+//                            -(i+1) 2
+//                      (q + 2      )  <= y.                     (2)
+//                        i
+//                                                            -(i+1)
+//      If (2) is false, then q   = q ; otherwise q   = q  + 2      .
+//                             i+1   i             i+1   i
+//
+//      With some algebric manipulation, it is not difficult to see
+//      that (2) is equivalent to
+//                             -(i+1)
+//                      s  +  2       <= y                       (3)
+//                       i                i
+//
+//      The advantage of (3) is that s  and y  can be computed by
+//                                    i      i
+//      the following recurrence formula:
+//          if (3) is false
+//
+//          s     =  s  ,       y    = y   ;                     (4)
+//           i+1      i          i+1    i
+//
+//      otherwise,
+//                         -i                      -(i+1)
+//          s     =  s  + 2  ,  y    = y  -  s  - 2              (5)
+//           i+1      i          i+1    i     i
+//
+//      One may easily use induction to prove (4) and (5).
+//      Note. Since the left hand side of (3) contain only i+2 bits,
+//            it does not necessary to do a full (53-bit) comparison
+//            in (3).
+//   3. Final rounding
+//      After generating the 53 bits result, we compute one more bit.
+//      Together with the remainder, we can decide whether the
+//      result is exact, bigger than 1/2ulp, or less than 1/2ulp
+//      (it will never equal to 1/2ulp).
+//      The rounding mode can be detected by checking whether
+//      huge + tiny is equal to huge, and whether huge - tiny is
+//      equal to huge for some floating point number "huge" and "tiny".
+//
+//
+// Notes:  Rounding mode detection omitted.  The constants "mask", "shift",
+// and "bias" are found in src/pkg/math/bits.go
+
+// Sqrt returns the square root of x.
+//
+// Special cases are:
+//	Sqrt(+Inf) = +Inf
+//	Sqrt(0) = 0
+//	Sqrt(x < 0) = NaN
+//	Sqrt(NaN) = NaN
+func sqrtGo(x float64) float64 {
+	// special cases
+	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
+	// when compiler does it for us
+	switch {
+	case x != x || x > MaxFloat64: // IsNaN(x) || IsInf(x, 1):
+		return x
+	case x == 0:
+		return 0
+	case x < 0:
+		return NaN()
+	}
+	ix := Float64bits(x)
+	// normalize x
+	exp := int((ix >> shift) & mask)
+	if exp == 0 { // subnormal x
+		for ix&1<<shift == 0 {
+			ix <<= 1
+			exp--
+		}
+		exp++
+	}
+	exp -= bias + 1 // unbias exponent
+	ix &^= mask << shift
+	ix |= 1 << shift
+	if exp&1 == 1 { // odd exp, double x to make it even
+		ix <<= 1
+	}
+	exp >>= 1 // exp = exp/2, exponent of square root
+	// generate sqrt(x) bit by bit
+	ix <<= 1
+	var q, s uint64               // q = sqrt(x)
+	r := uint64(1 << (shift + 1)) // r = moving bit from MSB to LSB
+	for r != 0 {
+		t := s + r
+		if t <= ix {
+			s = t + r
+			ix -= t
+			q += r
+		}
+		ix <<= 1
+		r >>= 1
+	}
+	// final rounding
+	if ix != 0 { // remainder, result not exact
+		q += q & 1 // round according to extra bit
+	}
+	ix = q>>1 + uint64(exp+bias)<<shift // significand + biased exponent
+	return Float64frombits(ix)
+}