#include <stdio.h>
#include <math.h>

static double my_fma(double x, double y, double z)
{
	return (long double) x*y + z;
}

int main(void)
{
	double x = 0x1p+1023, y = 0x1p+1, z = -0x1.fffffffffffffp+1023;

	printf("Computing\n"
	       "  %a * %a + %a\n"
	       "by double-precision fused muliply-add.\n\n"
	       "The result should be finite even though %a*%a overflows.\n\n",
		x, y, z, x, y);

	printf("unfused (this one should be infinite): %a\n", x*y+z);
	printf("GCC builtin: %a\n", __builtin_fma(x, y, z));
	printf("C library:   %a\n", (fma)(x, y, z));
	printf("Long double: %a\n", my_fma(x, y, z));

	return 0;
}