initial

1 files changed, 38 insertions, 0 deletions

diff --git a/libk/src/math/fma.c b/libk/src/math/fma.c
new file mode 100644
index 0000000..b26229e
--- /dev/null
+++ b/libk/src/math/fma.c
@@ -0,0 +1,38 @@
+#include <stdint.h>
+#include <math.h>
+
+double fma(double x, double y, double z) {
+	double xy, result;
+	union {double f; uint64_t i;} u;
+	int e;
+
+	xy = (double)x * y;
+	result = xy + z;
+	u.f = result;
+	e = u.i>>52 & 0x7ff;
+	if (
+        (u.i & 0x1fffffff) != 0x10000000 ||
+		e == 0x7ff ||
+		(result - xy == z && result - z == xy)
+    ) {
+		z = result;
+		return z;
+	}
+
+	double err;
+	int neg = u.i >> 63;
+	if (neg == (z > xy))
+		err = xy - result + z;
+	else
+		err = z - result + xy;
+	if (neg == (err < 0))
+		u.i++;
+	else
+		u.i--;
+	z = u.f;
+	return z;
+}
+
+float fmaf(float x, float y, float z) {
+    return fma(x, y, z);
+}