libk/src/math/fma.c


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
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);
}