| // 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 |
| |
| /* |
| Floating-point square root. |
| */ |
| |
| // 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 == 0 || x != x || x > MaxFloat64: // 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 + 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) |
| } |