golang sqrt 代码
文件路径：/src/math/sqrt.go
/ 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 algebraic 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 is 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/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 Sqrt(x float64) float64 {
if haveArchSqrt {
return archSqrt(x)
}
return sqrt(x)
}
/ Note: Sqrt is implemented in assembly on some systems.
/ Others have assembly stubs that jump to func sqrt below.
/ On systems where Sqrt is a single instruction, the compiler
/ may turn a direct call into a direct use of that instruction instead.
func sqrt(x float64) float64 {
/ special cases
switch {
case x == 0 || IsNaN(x) || IsInf(x, 1):
return x
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 / 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-1+bias)<<shift / significand + biased exponent
return Float64frombits(ix)
}