f8d9fa9e80
This upgrades all of libgo other than the runtime package to the Go 1.4 release. In Go 1.4 much of the runtime was rewritten into Go. Merging that code will take more time and will not change the API, so I'm putting it off for now. There are a few runtime changes anyhow, to accomodate other packages that rely on minor modifications to the runtime support. The compiler changes slightly to add a one-bit flag to each type descriptor kind that is stored directly in an interface, which for gccgo is currently only pointer types. Another one-bit flag (gcprog) is reserved because it is used by the gc compiler, but gccgo does not currently use it. There is another error check in the compiler since I ran across it during testing. gotools/: * Makefile.am (go_cmd_go_files): Sort entries. Add generate.go. * Makefile.in: Rebuild. From-SVN: r219627
149 lines
4.7 KiB
Go
149 lines
4.7 KiB
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
|
|
|
|
//extern sqrt
|
|
func libc_sqrt(float64) float64
|
|
|
|
func Sqrt(x float64) float64 {
|
|
return libc_sqrt(x)
|
|
}
|
|
|
|
// 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. 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 {
|
|
// 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)
|
|
}
|
|
|
|
func sqrtC(f float64, r *float64) {
|
|
*r = sqrt(f)
|
|
}
|