---
 mingw-w64-crt/Makefile.am |  4 ++--
 mingw-w64-crt/math/fma.S  | 42 --------------------------------
 mingw-w64-crt/math/fma.c  | 12 ++++++++++
 mingw-w64-crt/math/fmaf.S | 43 ---------------------------------
 mingw-w64-crt/math/fmaf.c | 10 ++++++++
 mingw-w64-crt/math/fmal.c | 61 ++++++++++++++++++++++++++++++++++++++++++++++-
 6 files changed, 84 insertions(+), 88 deletions(-)
 delete mode 100644 mingw-w64-crt/math/fma.S
 create mode 100644 mingw-w64-crt/math/fma.c
 delete mode 100644 mingw-w64-crt/math/fmaf.S
 create mode 100644 mingw-w64-crt/math/fmaf.c

diff --git a/mingw-w64-crt/Makefile.am b/mingw-w64-crt/Makefile.am
index 886fcf0..1c6e534 100644
--- a/mingw-w64-crt/Makefile.am
+++ b/mingw-w64-crt/Makefile.am
@@ -244,7 +244,6 @@ src_libmingwex=\
   \
   math/_chgsignl.S      math/ceil.S        math/ceilf.S          math/ceill.S  
       math/copysignl.S     \
   math/floor.S          math/floorf.S      math/floorl.S         \
-  math/fma.S            math/fmaf.S        \
   math/nearbyint.S      math/nearbyintf.S  math/nearbyintl.S     \
   math/trunc.S          math/truncf.S      \
   math/cbrt.c           \
@@ -252,7 +251,8 @@ src_libmingwex=\
   math/coshf.c          math/coshl.c       math/erfl.c           \
   math/expf.c           \
   math/fabs.c           math/fabsf.c       math/fabsl.c          math/fdim.c   
       math/fdimf.c         math/fdiml.c         \
-  math/fmal.c           math/fmax.c        math/fmaxf.c          math/fmaxl.c  
       math/fmin.c          math/fminf.c         \
+  math/fma.c            math/fmaf.c        math/fmal.c           \
+  math/fmax.c           math/fmaxf.c       math/fmaxl.c          math/fmin.c   
       math/fminf.c         \
   math/fminl.c          math/fp_consts.c   math/fp_constsf.c     \
   math/fp_constsl.c     math/fpclassify.c  math/fpclassifyf.c    
math/fpclassifyl.c   math/frexpf.c        \
   math/hypotf.c         math/hypot.c  math/hypotl.c      math/isnan.c          
math/isnanf.c        math/isnanl.c        \
diff --git a/mingw-w64-crt/math/fma.S b/mingw-w64-crt/math/fma.S
deleted file mode 100644
index 74becde..0000000
--- a/mingw-w64-crt/math/fma.S
+++ /dev/null
@@ -1,42 +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 <_mingw_mac.h>
-
-       .file   "fma.S"
-       .text
-#ifdef __x86_64__
-       .align 8
-#else
-       .align 4
-#endif
-       .p2align 4,,15
-       .globl __MINGW_USYMBOL(fma)
-       .def    __MINGW_USYMBOL(fma);   .scl    2;      .type   32;     .endef
-__MINGW_USYMBOL(fma):
-#if defined(_AMD64_) || defined(__x86_64__)
-       subq    $56, %rsp
-       movsd   %xmm0,(%rsp)
-       movsd   %xmm1,16(%rsp)
-       movsd   %xmm2,32(%rsp)
-       fldl    (%rsp)
-       fmull   16(%rsp)
-       fldl    32(%rsp)
-       faddp
-       fstpl   (%rsp)
-       movsd   (%rsp),%xmm0
-       addq    $56, %rsp
-       ret
-#elif defined(_ARM_) || defined(__arm__)
-       fmacd d2, d0, d1
-       fcpyd d0, d2
-       bx      lr
-#elif defined(_X86_) || defined(__i386__)
-       fldl    4(%esp)
-       fmull   12(%esp)
-       fldl    20(%esp)
-       faddp
-       ret
-#endif
diff --git a/mingw-w64-crt/math/fma.c b/mingw-w64-crt/math/fma.c
new file mode 100644
index 0000000..3703e00
--- /dev/null
+++ b/mingw-w64-crt/math/fma.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.
+ */
+long double fmal ( long double _x,  long double _y,  long double _z);
+
+double
+fma ( double _x,  double _y,  double _z)
+{
+  return (double)fmal(_x, _y, _z);
+}
diff --git a/mingw-w64-crt/math/fmaf.S b/mingw-w64-crt/math/fmaf.S
deleted file mode 100644
index 6bc7ef0..0000000
--- a/mingw-w64-crt/math/fmaf.S
+++ /dev/null
@@ -1,43 +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 <_mingw_mac.h>
-
-       .file   "fmaf.S"
-       .text
-#ifdef __x86_64__
-       .align 8
-#else
-       .align 2
-#endif
-       .p2align 4,,15
-       .globl __MINGW_USYMBOL(fmaf)
-       .def    __MINGW_USYMBOL(fmaf);  .scl    2;      .type   32;     .endef
-__MINGW_USYMBOL(fmaf):
-#if defined(_AMD64_) || defined(__x86_64__)
-       subq    $56, %rsp
-       movss   %xmm0,(%rsp)
-       movss   %xmm1,16(%rsp)
-       movss   %xmm2,32(%rsp)
-       flds    (%rsp)
-       fmuls   16(%rsp)
-       flds    32(%rsp)
-       faddp
-       fstps   (%rsp)
-       movss   (%rsp),%xmm0
-       addq    $56, %rsp
-       ret
-#elif defined(_ARM_) || defined(__arm__)
-       fmacs s2, s0, s1
-       fcpys s0, s2
-       bx      lr
-#elif defined(_X86_) || defined(__i386__)
-       flds    4(%esp)
-       fmuls   8(%esp)
-       flds    12(%esp)
-       faddp
-       ret
-#endif
-
diff --git a/mingw-w64-crt/math/fmaf.c b/mingw-w64-crt/math/fmaf.c
new file mode 100644
index 0000000..70a396e
--- /dev/null
+++ b/mingw-w64-crt/math/fmaf.c
@@ -0,0 +1,10 @@
+/**
+ * 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.
+ */
+float
+fmaf ( float _x,  float _y,  float _z)
+{
+  return (double)_x * (double)_y + (double)_z;
+}
diff --git a/mingw-w64-crt/math/fmal.c b/mingw-w64-crt/math/fmal.c
index 3d0eb02..5cf21f8 100644
--- a/mingw-w64-crt/math/fmal.c
+++ b/mingw-w64-crt/math/fmal.c
@@ -3,10 +3,69 @@
  * This file is part of the mingw-w64 runtime package.
  * No warranty is given; refer to the file DISCLAIMER.PD within this package.
  */
+#include <stdint.h>
+#include <stdbool.h>
+
 long double fmal ( long double _x,  long double _y,  long double _z);
 
+// 
https://en.wikipedia.org/wiki/Extended_precision#x86_Extended_Precision_Format
+typedef union x87reg_ {
+  long double f;
+  struct __attribute__((__packed__)) {
+    uint64_t f64 : 64; // fraction
+    uint16_t exp : 15; // exponent
+    uint16_t sgn :  1; // sign
+  };
+  struct __attribute__((__packed__)) {
+    uint32_t flo : 32;
+    uint32_t fhi : 32;
+//    uint16_t exp : 15;
+//    uint16_t sgn :  1;
+  };
+} x87reg;
+
+static inline void break_down(x87reg *restrict lo, x87reg *restrict hi, long 
double x){
+  hi->f = x;
+  const uint32_t flo = hi->flo;
+  const uint32_t exp = hi->exp;
+  const bool     sgn = hi->sgn;
+  // Erase low-order significant bits. `hi->f` now has only 32 significant 
bits.
+  hi->flo = 0;
+
+  if(flo == 0){
+    // If the low-order significant bits are all zeroes, return zero in 
`lo->f`.
+    lo->f64 = 0;
+    lo->exp = 0;
+  } else {
+    // How many bits should we shift to normalize the floating point value?
+    // We assume `long` has 32 bits here.
+    const unsigned shn = (unsigned)__builtin_clzl(flo) + 32;
+    if(shn < exp){
+      // `x` can be normalized!
+      lo->f64 = (uint64_t)flo << shn; // The MSB becomes one.
+      lo->exp = (exp - shn) & 0x7FFF; // Don'e forget to adjust the exponent.
+    } else {
+      // There are so few bits that `x` can't be normalized. Go with a 
denormal number.
+      lo->f64 = (uint64_t)flo << exp; // The MSB is zero.
+      lo->exp = 0;
+    }
+  }
+  lo->sgn = sgn;
+}
+static inline long double fpu_fma(long double x, long double y, long double z){
+  x87reg xlo, xhi, ylo, yhi;
+  break_down(&xlo, &xhi, x);
+  break_down(&ylo, &yhi, y);
+  // The order of these four statements is essential. Don't move them around.
+  long double ret = z;
+  ret += xhi.f * yhi.f;                 // This is the most significant item.
+  ret += xhi.f * ylo.f + xlo.f * yhi.f; // They are equally significant.
+  ret += xlo.f * ylo.f;                 // This is the least significant item.
+  return ret;
+}
+
 long double
 fmal ( long double _x,  long double _y,  long double _z)
 {
-  return ((_x * _y) + _z);
+  return fpu_fma(_x, _y, _z);
 }
-- 
2.9.1

------------------------------------------------------------------------------
_______________________________________________
Mingw-w64-public mailing list
Mingw-w64-public@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public

Reply via email to