#include <stdio.h>
#include <float.h>
#include <math.h>
#include <fenv.h>
#include <errno.h>

double
pg_hypot(double x, double y)
{
	double		yx;

	/* Handle INF and NaN properly */
	if (isinf(x) || isinf(y))
		return (double) INFINITY;

	if (isnan(x) || isnan(y))
		return (double) NAN;

	/* Else, drop any minus signs */
	x = fabs(x);
	y = fabs(y);

	/* Swap x and y if needed to make x the larger one */
	if (x < y)
	{
		double		temp = x;

		x = y;
		y = temp;
	}

	/*
	 * If y is zero, the hypotenuse is x.  This test saves a few cycles in
	 * such cases, but more importantly it also protects against
	 * divide-by-zero errors, since now x >= y.
	 */
	if (y == 0.0)
		return x;

	/* Determine the hypotenuse */
	yx = y / x;
	return x * sqrt(1.0 + (yx * yx));
}

void
setfeflags(int *invalid, int *divzero, int *overflow, int *underflow)
{
	int err = fetestexcept(FE_INVALID | FE_DIVBYZERO |
						   FE_OVERFLOW | FE_UNDERFLOW);
	*invalid = ((err & FE_INVALID) != 0);
	*divzero = ((err & FE_DIVBYZERO) != 0);
	*overflow = ((err & FE_OVERFLOW) != 0);
	*underflow = ((err & FE_UNDERFLOW) != 0);
}

int
main(void)
{
	double x;
	double y;
	double p[][2] =
		{
			{2.2e-308, 2.2e-308},
			{2.2e-308, 1.7e307},
			{1.7e307, 1.7e307},
			{1.7e308, 1.7e308},
			{2.2e-308, DBL_MAX},
			{1.7e308, DBL_MAX},
			{DBL_MIN, DBL_MAX},
			{DBL_MAX, DBL_MAX},
			{1.7e307, INFINITY},
			{2.2e-308, INFINITY},
			{0, INFINITY},
			{DBL_MIN, INFINITY},
			{INFINITY, INFINITY},
			{1, NAN},
			{INFINITY, NAN},
			{NAN, NAN},
			{0.0, 0.0}
		};
	int i;


	for (i = 0 ; p[i][0] != 0.0 || p[i][1] != 0.0 ; i++)
	{
		int errno1, errno2;
		int invalid1 = 0, divzero1 = 0, overflow1 = 0, underflow1 = 0;
		int invalid2 = 0, divzero2 = 0, overflow2 = 0, underflow2 = 0;
		int cmp;

		errno = 0;
		feclearexcept(FE_ALL_EXCEPT);
		x = hypot(p[i][0], p[i][1]);
		errno1 = errno;
		setfeflags(&invalid1, &divzero1, &overflow1, &underflow1);

		errno = 0;
		feclearexcept(FE_ALL_EXCEPT);
		y = pg_hypot(p[i][0], p[i][1]);
		errno2 = errno;
		setfeflags(&invalid2, &divzero2, &overflow2, &underflow2);

		/* assume NaN == NaN here */
		if (isnan(x) && isnan(y))
			cmp = 1;
		else
			cmp = (x == y);

		printf("%d: hypot=%.16lg (== 0.0 is %d, < DBL_MIN is %d)\n  pg_hypot=%.16lg (== 0.0 is %d, < DBL_MIN is %d)\n  equal=%d, hypot(errno:%d, inval:%d, div0:%d, of=%d, uf=%d), pg_hypot(errno:%d, inval:%d, div0:%d, of=%d, uf=%d)\n",
			   i, x, x == 0.0, x < DBL_MIN, y, y == 0.0, y < DBL_MIN, cmp,
			   errno1, invalid1, divzero1, overflow1, underflow1,
			   errno2, invalid2, divzero2, overflow2, underflow2);
	}
	return 0;
}