crt: arm: Use a bit twiddling implementation of rint*
Split the implementation file per architecture, and use a BSD/Sun
implementation of rint/rintf for arm, as the previous arm
implementation used signed 32 bit integer intermediates.
Signed-off-by: Martin Storsjö <martin@martin.st>
diff --git a/mingw-w64-crt/math/arm/s_rint.c b/mingw-w64-crt/math/arm/s_rint.c
new file mode 100644
index 0000000..823a2e0
--- /dev/null
+++ b/mingw-w64-crt/math/arm/s_rint.c
@@ -0,0 +1,86 @@
+/* @(#)s_rint.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+
+/*
+ * rint(x)
+ * Return x rounded to integral value according to the prevailing
+ * rounding mode.
+ * Method:
+ * Using floating addition.
+ * Exception:
+ * Inexact flag raised if x not equal to rint(x).
+ */
+
+#include <float.h>
+
+#include "../bsd_private_base.h"
+
+static const double
+TWO52[2]={
+ 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
+ -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
+};
+
+double
+rint(double x)
+{
+ int32_t i0,j0,sx;
+ u_int32_t i,i1;
+ double w,t;
+ EXTRACT_WORDS(i0,i1,x);
+ sx = (i0>>31)&1;
+ j0 = ((i0>>20)&0x7ff)-0x3ff;
+ if(j0<20) {
+ if(j0<0) {
+ if(((i0&0x7fffffff)|i1)==0) return x;
+ i1 |= (i0&0x0fffff);
+ i0 &= 0xfffe0000;
+ i0 |= ((i1|-i1)>>12)&0x80000;
+ SET_HIGH_WORD(x,i0);
+ STRICT_ASSIGN(double,w,TWO52[sx]+x);
+ t = w-TWO52[sx];
+ GET_HIGH_WORD(i0,t);
+ SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
+ return t;
+ } else {
+ i = (0x000fffff)>>j0;
+ if(((i0&i)|i1)==0) return x; /* x is integral */
+ i>>=1;
+ if(((i0&i)|i1)!=0) {
+ /*
+ * Some bit is set after the 0.5 bit. To avoid the
+ * possibility of errors from double rounding in
+ * w = TWO52[sx]+x, adjust the 0.25 bit to a lower
+ * guard bit. We do this for all j0<=51. The
+ * adjustment is trickiest for j0==18 and j0==19
+ * since then it spans the word boundary.
+ */
+ if(j0==19) i1 = 0x40000000; else
+ if(j0==18) i1 = 0x80000000; else
+ i0 = (i0&(~i))|((0x20000)>>j0);
+ }
+ }
+ } else if (j0>51) {
+ if(j0==0x400) return x+x; /* inf or NaN */
+ else return x; /* x is integral */
+ } else {
+ i = ((u_int32_t)(0xffffffff))>>(j0-20);
+ if((i1&i)==0) return x; /* x is integral */
+ i>>=1;
+ if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
+ }
+ INSERT_WORDS(x,i0,i1);
+ STRICT_ASSIGN(double,w,TWO52[sx]+x);
+ return w-TWO52[sx];
+}
diff --git a/mingw-w64-crt/math/arm/s_rintf.c b/mingw-w64-crt/math/arm/s_rintf.c
new file mode 100644
index 0000000..0ec0170
--- /dev/null
+++ b/mingw-w64-crt/math/arm/s_rintf.c
@@ -0,0 +1,51 @@
+/* s_rintf.c -- float version of s_rint.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+
+#include <float.h>
+#include <stdint.h>
+
+#include "../bsd_private_base.h"
+
+static const float
+TWO23[2]={
+ 8.3886080000e+06, /* 0x4b000000 */
+ -8.3886080000e+06, /* 0xcb000000 */
+};
+
+float
+rintf(float x)
+{
+ int32_t i0,j0,sx;
+ float w,t;
+ GET_FLOAT_WORD(i0,x);
+ sx = (i0>>31)&1;
+ j0 = ((i0>>23)&0xff)-0x7f;
+ if(j0<23) {
+ if(j0<0) {
+ if((i0&0x7fffffff)==0) return x;
+ STRICT_ASSIGN(float,w,TWO23[sx]+x);
+ t = w-TWO23[sx];
+ GET_FLOAT_WORD(i0,t);
+ SET_FLOAT_WORD(t,(i0&0x7fffffff)|(sx<<31));
+ return t;
+ }
+ STRICT_ASSIGN(float,w,TWO23[sx]+x);
+ return w-TWO23[sx];
+ }
+ if(j0==0x80) return x+x; /* inf or NaN */
+ else return x; /* x is integral */
+}
diff --git a/mingw-w64-crt/math/arm64/rint.c b/mingw-w64-crt/math/arm64/rint.c
new file mode 100644
index 0000000..03a7ffa
--- /dev/null
+++ b/mingw-w64-crt/math/arm64/rint.c
@@ -0,0 +1,12 @@
+/**
+ * This file has no copyright assigned and is placed in the Public Domain.
+ * This file is part of the mingw-w64 runtime package.
+ * No warranty is given; refer to the file DISCLAIMER.PD within this package.
+ */
+#include <math.h>
+
+double rint (double x) {
+ double retval = 0.0;
+ __asm__ __volatile__ ("frintx %d0, %d1\n\t" : "=w" (retval) : "w" (x));
+ return retval;
+}
diff --git a/mingw-w64-crt/math/arm64/rintf.c b/mingw-w64-crt/math/arm64/rintf.c
new file mode 100644
index 0000000..914769a
--- /dev/null
+++ b/mingw-w64-crt/math/arm64/rintf.c
@@ -0,0 +1,12 @@
+/**
+ * This file has no copyright assigned and is placed in the Public Domain.
+ * This file is part of the mingw-w64 runtime package.
+ * No warranty is given; refer to the file DISCLAIMER.PD within this package.
+ */
+#include <math.h>
+
+float rintf (float x) {
+ float retval = 0.0F;
+ __asm__ __volatile__ ("frintx %s0, %s1\n\t" : "=w" (retval) : "w" (x));
+ return retval;
+}
diff --git a/mingw-w64-crt/math/bsd_private_base.h b/mingw-w64-crt/math/bsd_private_base.h
index f31c75d..de5b935 100644
--- a/mingw-w64-crt/math/bsd_private_base.h
+++ b/mingw-w64-crt/math/bsd_private_base.h
@@ -10,6 +10,7 @@
*/
#include <inttypes.h>
+#include <float.h>
typedef unsigned int u_int32_t;
@@ -99,3 +100,24 @@
gf_u.word = (i); \
(d) = gf_u.value; \
} while(0)
+
+
+#ifdef FLT_EVAL_METHOD
+/*
+ * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
+ */
+#if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
+#define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval))
+#else
+#define STRICT_ASSIGN(type, lval, rval) do { \
+ volatile type __lval; \
+ \
+ if (sizeof(type) >= sizeof(long double)) \
+ (lval) = (rval); \
+ else { \
+ __lval = (rval); \
+ (lval) = __lval; \
+ } \
+} while (0)
+#endif
+#endif /* FLT_EVAL_METHOD */
diff --git a/mingw-w64-crt/math/rint.c b/mingw-w64-crt/math/rint.c
deleted file mode 100644
index d883314..0000000
--- a/mingw-w64-crt/math/rint.c
+++ /dev/null
@@ -1,24 +0,0 @@
-/**
- * This file has no copyright assigned and is placed in the Public Domain.
- * This file is part of the mingw-w64 runtime package.
- * No warranty is given; refer to the file DISCLAIMER.PD within this package.
- */
-#include <math.h>
-
-double rint (double x) {
- double retval = 0.0;
-#if defined(_AMD64_) || defined(__x86_64__) || defined(_X86_) || defined(__i386__)
- __asm__ __volatile__ ("frndint;" : "=t" (retval) : "0" (x));
-#elif defined(__arm__) || defined(_ARM_)
- if (isnan(x) || isinf(x))
- return x;
- float temp;
- __asm__ __volatile__ (
- "vcvtr.s32.f64 %[tmp], %[src]\n\t"
- "vcvt.f64.s32 %[dst], %[tmp]\n\t"
- : [dst] "=w" (retval), [tmp] "=t" (temp) : [src] "w" (x));
-#elif defined(__aarch64__) || defined(_ARM64_)
- __asm__ __volatile__ ("frintx %d0, %d1\n\t" : "=w" (retval) : "w" (x));
-#endif
- return retval;
-}
diff --git a/mingw-w64-crt/math/rintf.c b/mingw-w64-crt/math/rintf.c
deleted file mode 100644
index bfae515..0000000
--- a/mingw-w64-crt/math/rintf.c
+++ /dev/null
@@ -1,23 +0,0 @@
-/**
- * This file has no copyright assigned and is placed in the Public Domain.
- * This file is part of the mingw-w64 runtime package.
- * No warranty is given; refer to the file DISCLAIMER.PD within this package.
- */
-#include <math.h>
-
-float rintf (float x) {
- float retval = 0.0F;
-#if defined(_AMD64_) || defined(__x86_64__) || defined(_X86_) || defined(__i386__)
- __asm__ __volatile__ ("frndint;": "=t" (retval) : "0" (x));
-#elif defined(__arm__) || defined(_ARM_)
- if (isnan(x) || isinf(x))
- return x;
- __asm__ __volatile__ (
- "vcvtr.s32.f32 %[dst], %[src]\n\t"
- "vcvt.f32.s32 %[dst], %[dst]\n\t"
- : [dst] "=t" (retval) : [src] "w" (x));
-#elif defined(__aarch64__) || defined(_ARM64_)
- __asm__ __volatile__ ("frintx %s0, %s1\n\t" : "=w" (retval) : "w" (x));
-#endif
- return retval;
-}
diff --git a/mingw-w64-crt/math/x86/rint.c b/mingw-w64-crt/math/x86/rint.c
new file mode 100644
index 0000000..55a73d1
--- /dev/null
+++ b/mingw-w64-crt/math/x86/rint.c
@@ -0,0 +1,12 @@
+/**
+ * This file has no copyright assigned and is placed in the Public Domain.
+ * This file is part of the mingw-w64 runtime package.
+ * No warranty is given; refer to the file DISCLAIMER.PD within this package.
+ */
+#include <math.h>
+
+double rint (double x) {
+ double retval = 0.0;
+ __asm__ __volatile__ ("frndint;" : "=t" (retval) : "0" (x));
+ return retval;
+}
diff --git a/mingw-w64-crt/math/x86/rintf.c b/mingw-w64-crt/math/x86/rintf.c
new file mode 100644
index 0000000..f5a6bf2
--- /dev/null
+++ b/mingw-w64-crt/math/x86/rintf.c
@@ -0,0 +1,12 @@
+/**
+ * This file has no copyright assigned and is placed in the Public Domain.
+ * This file is part of the mingw-w64 runtime package.
+ * No warranty is given; refer to the file DISCLAIMER.PD within this package.
+ */
+#include <math.h>
+
+float rintf (float x) {
+ float retval = 0.0F;
+ __asm__ __volatile__ ("frndint;": "=t" (retval) : "0" (x));
+ return retval;
+}
