Note: sqrt_decl.go already in src/pkg/math/.
R=rsc
CC=golang-dev
https://golang.org/cl/183155
OFILES_amd64=\
sqrt_amd64.$O\
+OFILES_386=\
+ sqrt_386.$O\
+
OFILES=\
$(OFILES_$(GOARCH))
sin.go\
sinh.go\
sqrt.go\
+ sqrt_port.go\
tan.go\
tanh.go\
unsafe.go\
Acos(.5)
}
}
+
+func BenchmarkSqrt(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Sqrt(10)
+ }
+}
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(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) }
--- /dev/null
+// 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
--- /dev/null
+// 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)
+}