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 */
6 float fmaf(float x, float y, float z);
7
8 #if defined(_ARM_) || defined(__arm__)
9
10 /* Use hardware FMA on ARM. */
fmaf(float x,float y,float z)11 float fmaf(float x, float y, float z){
12 __asm__ (
13 "fmacs %0, %1, %2 \n"
14 : "+t"(z)
15 : "t"(x), "t"(y)
16 );
17 return z;
18 }
19
20 #elif defined(_ARM64_) || defined(__aarch64__)
21
22 /* Use hardware FMA on ARM64. */
fmaf(float x,float y,float z)23 float fmaf(float x, float y, float z){
24 __asm__ (
25 "fmadd %s0, %s1, %s2, %s0 \n"
26 : "+w"(z)
27 : "w"(x), "w"(y)
28 );
29 return z;
30 }
31
32 #elif defined(_AMD64_) || defined(__x86_64__) || defined(_X86_) || defined(__i386__)
33
34 #include <math.h>
35 #include <stdint.h>
36
37 /* This is in accordance with the IEC 559 single-precision format.
38 * Be advised that due to the hidden bit, the higher half actually has 11 bits.
39 * Multiplying two 13-bit numbers will cause a 1-ULP error, which we cannot
40 * avoid. It is kept in the very last position.
41 */
42 typedef union iec559_float_ {
43 struct __attribute__((__packed__)) {
44 uint32_t mlo : 13;
45 uint32_t mhi : 10;
46 uint32_t exp : 8;
47 uint32_t sgn : 1;
48 };
49 float f;
50 } iec559_float;
51
break_down(iec559_float * restrict lo,iec559_float * restrict hi,float x)52 static inline void break_down(iec559_float *restrict lo, iec559_float *restrict hi, float x) {
53 hi->f = x;
54 /* Erase low-order significant bits. `hi->f` now has only 11 significant bits. */
55 hi->mlo = 0;
56 /* Store the low-order half. It will be normalized by the hardware. */
57 lo->f = x - hi->f;
58 /* Preserve signness in case of zero. */
59 lo->sgn = hi->sgn;
60 }
61
fmaf(float x,float y,float z)62 float fmaf(float x, float y, float z) {
63 /*
64 POSIX-2013:
65 1. If x or y are NaN, a NaN shall be returned.
66 2. If x multiplied by y is an exact infinity and z is also an infinity
67 but with the opposite sign, a domain error shall occur, and either a NaN
68 (if supported), or an implementation-defined value shall be returned.
69 3. If one of x and y is infinite, the other is zero, and z is not a NaN,
70 a domain error shall occur, and either a NaN (if supported), or an
71 implementation-defined value shall be returned.
72 4. If one of x and y is infinite, the other is zero, and z is a NaN, a NaN
73 shall be returned and a domain error may occur.
74 5. If x* y is not 0*Inf nor Inf*0 and z is a NaN, a NaN shall be returned.
75 */
76 /* Check whether the result is finite. */
77 float ret = x * y + z;
78 if(!isfinite(ret)) {
79 return ret; /* If this naive check doesn't yield a finite value, the FMA isn't
80 likely to return one either. Forward the value as is. */
81 }
82 iec559_float xlo, xhi, ylo, yhi;
83 break_down(&xlo, &xhi, x);
84 break_down(&ylo, &yhi, y);
85 /* The order of these four statements is essential. Don't move them around. */
86 ret = z;
87 ret += xhi.f * yhi.f; /* The most significant item comes first. */
88 ret += xhi.f * ylo.f + xlo.f * yhi.f; /* They are equally significant. */
89 ret += xlo.f * ylo.f; /* The least significant item comes last. */
90 return ret;
91 }
92
93 #else
94
95 #error Please add FMA implementation for this platform.
96
97 #endif
98