From ee7bb07a53da1c400f4e1130517c362e302be212 Mon Sep 17 00:00:00 2001 From: Russ Cox Date: Mon, 12 May 2014 10:55:33 -0400 Subject: [PATCH] runtime: add copy of math.sqrt for use by arm softfloat If it's not used (such as on other systems or if softfloat is disabled) the linker will discard it. The alternative is to teach cmd/go that every binary depends on math implicitly on arm. I started down that path but it's too scary. If we're going to get dependencies right we should get dependencies right. Fixes #6994. LGTM=bradfitz, dave R=golang-codereviews, bradfitz, dave CC=golang-codereviews https://golang.org/cl/95290043 --- src/pkg/runtime/softfloat_arm.c | 4 +- src/pkg/runtime/sqrt.go | 150 ++++++++++++++++++++++++++++++++ 2 files changed, 152 insertions(+), 2 deletions(-) create mode 100644 src/pkg/runtime/sqrt.go diff --git a/src/pkg/runtime/softfloat_arm.c b/src/pkg/runtime/softfloat_arm.c index f5801dde43..29a52bd0e4 100644 --- a/src/pkg/runtime/softfloat_arm.c +++ b/src/pkg/runtime/softfloat_arm.c @@ -16,7 +16,7 @@ #define FLAGS_V (1U << 28) void runtime·abort(void); -void math·sqrtC(uint64, uint64*); +void runtime·sqrtC(uint64, uint64*); static uint32 trace = 0; @@ -413,7 +413,7 @@ stage3: // regd, regm are 4bit variables break; case 0xeeb10bc0: // D[regd] = sqrt D[regm] - math·sqrtC(getd(regm), &uval); + runtime·sqrtC(getd(regm), &uval); putd(regd, uval); if(trace) diff --git a/src/pkg/runtime/sqrt.go b/src/pkg/runtime/sqrt.go new file mode 100644 index 0000000000..34a8c3806b --- /dev/null +++ b/src/pkg/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<>= 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)<