#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);
}