diff options
Diffstat (limited to 'src/runtime/sqrt.go')
| -rw-r--r-- | src/runtime/sqrt.go | 150 |
1 files changed, 150 insertions, 0 deletions
diff --git a/src/runtime/sqrt.go b/src/runtime/sqrt.go new file mode 100644 index 000000000..34a8c3806 --- /dev/null +++ b/src/runtime/sqrt.go @@ -0,0 +1,150 @@ +// 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) +} |