diff options
Diffstat (limited to 'src/pkg/runtime/sqrt.go')
| -rw-r--r-- | src/pkg/runtime/sqrt.go | 150 |
1 files changed, 0 insertions, 150 deletions
diff --git a/src/pkg/runtime/sqrt.go b/src/pkg/runtime/sqrt.go deleted file mode 100644 index 34a8c3806..000000000 --- a/src/pkg/runtime/sqrt.go +++ /dev/null @@ -1,150 +0,0 @@ -// 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. - -// Copy of math/sqrt.go, here for use by ARM softfloat. - -package runtime - -import "unsafe" - -// 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 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. - -const ( - uvnan = 0x7FF8000000000001 - uvinf = 0x7FF0000000000000 - uvneginf = 0xFFF0000000000000 - mask = 0x7FF - shift = 64 - 11 - 1 - bias = 1023 - maxFloat64 = 1.797693134862315708145274237317043567981e+308 // 2**1023 * (2**53 - 1) / 2**52 -) - -func float64bits(f float64) uint64 { return *(*uint64)(unsafe.Pointer(&f)) } -func float64frombits(b uint64) float64 { return *(*float64)(unsafe.Pointer(&b)) } - -func sqrt(x float64) float64 { - // special cases - switch { - case x == 0 || x != x || x > maxFloat64: - 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) -} - -func sqrtC(f float64, r *float64) { - *r = sqrt(f) -} |