qbit/go
1
0
Fork 0
mirror of https://github.com/golang/go synced 2024-11-25 01:57:56 -07:00
Code Releases Activity

math: Sqrt using 386 FPU.

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

R=rsc
CC=golang-dev
https://golang.org/cl/183155
This commit is contained in:
Charles L. Dorian 2010-01-10 15:41:07 -08:00 committed by Russ Cox
parent 5328df6534
commit 5336cd8f91
5 changed files with 162 additions and 129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/pkg/math/Makefile
View File

@ -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\

6
src/pkg/math/all_test.go
View File

@ -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)
}
}

130
src/pkg/math/sqrt.go
View File

@ -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) }

10
src/pkg/math/sqrt_386.s Normal file
View File

@ -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

141
src/pkg/math/sqrt_port.go Normal file
View File

@ -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)
}