1/**
 2 * This file has no copyright assigned and is placed in the Public Domain.
 3 * This file is part of the mingw-w64 runtime package.
 4 * No warranty is given; refer to the file DISCLAIMER.PD within this package.
 5 */
 6long double fmal(long double x, long double y, long double z);
 7
 8#if __SIZEOF_LONG_DOUBLE__ == __SIZEOF_DOUBLE__
 9
10double fma(double x, double y, double z);
11
12/* On ARM `long double` is 64 bits. And ARM has hardware FMA. */
13long double fmal(long double x, long double y, long double z){
14  return fma(x, y, z);
15}
16
17#elif defined(_AMD64_) || defined(__x86_64__) || defined(_X86_) || defined(__i386__)
18
19/**
20 * x87-specific software-emulated FMA by LH_Mouse (lh_mouse at 126 dot com).
21 * This file is donated to the mingw-w64 project.
22 * Note: This file requires C99 support to compile.
23 */
24
25#include <math.h>
26#include <stdint.h>
27
28/* See <https://en.wikipedia.org/wiki/Extended_precision#x86_extended_precision_format>.
29 * Note the higher half of the mantissa has fewer significant bits than the lower
30 * half, which reduces rounding errors in the more significant position but increases
31 * them in the other end.
32 */
33typedef union x87reg_ {
34  struct __attribute__((__packed__)) {
35    uint64_t mlo : 33;
36    uint64_t mhi : 31;
37    uint16_t exp : 15;
38    uint16_t sgn :  1;
39  };
40  long double f;
41} x87reg;
42
43static inline void break_down(x87reg *restrict lo, x87reg *restrict hi, long double x) {
44  hi->f = x;
45  /* Erase low-order significant bits. `hi->f` now has only 31 significant bits. */
46  hi->mlo = 0;
47  /* Store the low-order half. It will be normalized by the hardware. */
48  lo->f = x - hi->f;
49  /* Preserve signness in case of zero. */
50  lo->sgn = hi->sgn;
51}
52
53long double fmal(long double x, long double y, long double z) {
54  /*
55    POSIX-2013:
56    1. If x or y are NaN, a NaN shall be returned.
57    2. If x multiplied by y is an exact infinity and z is also an infinity
58       but with the opposite sign, a domain error shall occur, and either a NaN
59       (if supported), or an implementation-defined value shall be returned.
60    3. If one of x and y is infinite, the other is zero, and z is not a NaN,
61       a domain error shall occur, and either a NaN (if supported), or an
62       implementation-defined value shall be returned.
63    4. If one of x and y is infinite, the other is zero, and z is a NaN, a NaN
64       shall be returned and a domain error may occur.
65    5. If x* y is not 0*Inf nor Inf*0 and z is a NaN, a NaN shall be returned.
66  */
67  /* Check whether the result is finite. */
68  long double ret = x * y + z;
69  if(!isfinite(ret)) {
70    return ret; /* If this naive check doesn't yield a finite value, the FMA isn't
71                   likely to return one either. Forward the value as is. */
72  }
73  x87reg xlo, xhi, ylo, yhi;
74  break_down(&xlo, &xhi, x);
75  break_down(&ylo, &yhi, y);
76  /* The order of these four statements is essential. Don't move them around. */
77  ret = z;
78  ret += xhi.f * yhi.f;                 /* The most significant item comes first. */
79  ret += xhi.f * ylo.f + xlo.f * yhi.f; /* They are equally significant. */
80  ret += xlo.f * ylo.f;                 /* The least significant item comes last. */
81  return ret;
82}
83
84#else
85
86#error Please add FMA implementation for this platform.
87
88#endif