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/2119
Reviewed-by: David Crawshaw <crawshaw@golang.org>
Reviewed-by: Dave Cheney <dave@cheney.net>

diff --git a/src/lib9/math_darwin_arm.c b/src/lib9/math_darwin_arm.c
new file mode 100644 (file)
index 0000000..ff5b0ce
--- /dev/null
+++ b/src/lib9/math_darwin_arm.c
@@ -0,0 +1,102 @@
+// 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"
+                    );
+}