lib9: provide alterantive implementation of ldexp and frexp

Libc on Darwin/ARM has a buggy implementation of ldexp and frexp
that could not handle denormal numbers.

Also disable VFP runfast (flush-to-zero) mode so that the gc
compiler can correctly handle denormal constants used in math and
strconv tests.

Change-Id: Ie64220b882f414e0b37f406f38181c3586104d46
Reviewed-on: https://go-review.googlesource.com/2119Reviewed-by: David Crawshaw <crawshaw@golang.org>
Reviewed-by: Dave Cheney <dave@cheney.net>

src/lib9/math_darwin_arm.c

+// Copyright 2014 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.
+
+// Some Darwin/ARM libc versions fail to provide a standard compliant version
+// of frexp and ldexp that could handle denormal floating point numbers.
+// The frexp and ldexp implementations are translated from their Go version.
+
+#include <stdint.h>
+
+// Assume double and uint64_t are using the same endian.
+union dint64 {
+	double d;
+	uint64_t u;
+};
+
+static const uint64_t mask = 0x7FF, bias = 1023;
+static const int shift = 64 - 11 - 1;
+static const uint64_t uvnan = 0x7FF8000000000001ULL, uvinf = 0x7FF0000000000000ULL,
+	  uvneginf = 0xFFF0000000000000ULL;
+static const double smallestnormal = 2.2250738585072014e-308; // 2**-1022
+
+static inline uint64_t float64bits(double x) {
+	union dint64 u;
+	u.d = x;
+	return u.u;
+}
+static inline double float64frombits(uint64_t x) {
+	union dint64 u;
+	u.u = x;
+	return u.d;
+}
+static inline int isinf(double x) {
+	return float64bits(x) == uvinf || float64bits(x) == uvneginf;
+}
+static inline int isnan(double x) {
+	return x != x;
+}
+extern double fabs(double);
+static double normalize(double x, int *exp) {
+	if (fabs(x) < smallestnormal) {
+		*exp = -52;
+		return x * (double)(1LL<<52);
+	}
+	*exp = 0;
+	return x;
+}
+
+double ldexp(double frac, int exp) {
+	// special cases
+	if (frac == 0.0) return frac;
+	if (isinf(frac) || isnan(frac)) return frac;
+
+	int e;
+	frac = normalize(frac, &e);
+	exp += e;
+	uint64_t x = float64bits(frac);
+	exp += (int)((x>>shift)&mask) - bias;
+	if (exp < -1074) { // underflow
+		if (frac < 0.0) return float64frombits(1ULL<<63); // -0.0
+		return 0.0;
+	}
+	if (exp > 1023) { // overflow
+		if (frac < 0.0) return float64frombits(uvneginf);
+		return float64frombits(uvinf);
+	}
+	double m = 1;
+	if (exp < -1022) { // denormal
+		exp += 52;
+		m = 1.0 / (double)(1ULL<<52);
+	}
+	x &= ~(mask << shift);
+	x |= (uint64_t)(exp+bias) << shift;
+	return m * float64frombits(x);
+}
+
+double frexp(double f, int *exp) {
+	*exp = 0;
+	// special cases
+	if (f == 0.0) return f;
+	if (isinf(f) || isnan(f)) return f;
+
+	f = normalize(f, exp);
+	uint64_t x = float64bits(f);
+	*exp += (int)((x>>shift)&mask) - bias + 1;
+	x &= ~(mask << shift);
+	x |= (-1 + bias) << shift;
+	return float64frombits(x);
+}
+
+// On Darwin/ARM, the kernel insists on running VFP in runfast mode, and it
+// cannot deal with denormal floating point numbers in that mode, so we have
+// to disable the runfast mode if the client uses ldexp/frexp (i.e. 5g).
+void disable_vfp_runfast(void) __attribute__((constructor));
+void disable_vfp_runfast(void) {
+    __asm__ volatile (
+		      "fmrx r0, fpscr\n"
+		      "bic r0, r0, $0x03000000\n"
+		      "fmxr fpscr, r0\n"
+		      : : : "r0"
+		     );
+}