diff options
Diffstat (limited to 'src/pkg/math/sqrt.go')
| -rw-r--r-- | src/pkg/math/sqrt.go | 163 |
1 files changed, 120 insertions, 43 deletions
diff --git a/src/pkg/math/sqrt.go b/src/pkg/math/sqrt.go index 1e2209f2a..a3a3119fe 100644 --- a/src/pkg/math/sqrt.go +++ b/src/pkg/math/sqrt.go @@ -4,13 +4,83 @@ package math - -/* - * sqrt returns the square root of its floating - * point argument. Newton's method. - * - * calls frexp - */ +// 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. // @@ -18,48 +88,55 @@ package math // Sqrt(+Inf) = +Inf // Sqrt(0) = 0 // Sqrt(x < 0) = NaN +// Sqrt(NaN) = NaN func Sqrt(x float64) float64 { - if IsInf(x, 1) { + // 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 - } - - if x <= 0 { - if x < 0 { - return NaN() - } + case x == 0: return 0 + case x < 0: + return NaN() } - - y, exp := Frexp(x) - for y < 0.5 { - y = y * 2 - exp = exp - 1 - } - - if exp&1 != 0 { - y = y * 2 - exp = exp - 1 - } - temp := 0.5 * (1 + y) - - for exp > 60 { - temp = temp * float64(1<<30) - exp = exp - 60 + ix := Float64bits(x) + // normalize x + exp := int((ix >> shift) & mask) + if exp == 0 { // subnormal x + for ix&1<<shift == 0 { + ix <<= 1 + exp-- + } + exp++ } - for exp < -60 { - temp = temp / float64(1<<30) - exp = exp + 60 + 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 } - if exp >= 0 { - exp = 1 << uint(exp/2) - temp = temp * float64(exp) - } else { - exp = 1 << uint(-exp/2) - temp = temp / float64(exp) + 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 } - - for i := 0; i <= 4; i++ { - temp = 0.5 * (temp + x/temp) + // final rounding + if ix != 0 { // remainder, result not exact + q += q & 1 // round according to extra bit } - return temp + ix = q>>1 + 0x3fe0000000000000 // q/2 + 0.5 + ix += uint64(exp) << shift + return Float64frombits(ix) } |