/*
 * Verify the functioning of hypot(3) on integer results
 */

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

#define sqr(x) ((x) * (x))
#define my_hypot(x, y) (sqrt(sqr((x)) + sqr((y))))

#define RESOLUTION 100

double ulp(double);

#ifndef USER
main()
{
	int p, q;

	// Generate Pythagorean triplets; see http://www.friesian.com/pythag.htm
	for (p = 1; p < RESOLUTION; p += 2)
		for (q = 1; q < p; q += 2) {
			double x = p * q;
			double y = (sqr(p) - sqr(q)) / 2.;
			double h = (sqr(p) + sqr(q)) / 2.;
			if (my_hypot(x, y) != h || hypot(x, y) != h) {
				double ulps = fabs(hypot(x, y) - my_hypot(x, y)) / ulp(hypot(x, y));
#ifdef ASCII
				printf("x=%.99g y=%.99g h=%.99g ulp=%g hm=%.99g hr=%.99g\n", x, y, h, ulps, my_hypot(x,y), hypot(x, y));
#else
				printf("<tr><td>%g</td><td>%g</td><td>%.99g</td><td>%.99g</td><td>%g</td></tr>\n", x, y, my_hypot(x,y), hypot(x, y), ulps);
#endif
			}
	}
	return 0;
}
#endif

#ifdef USER
main(int argc, char *argv[])
{
	int i, j;

	double x = atoi(argv[1]);
	double y = atoi(argv[2]);
	double ulps = fabs(hypot(x, y) - my_hypot(x, y)) / ulp(hypot(x, y));
	if (ulps > 0)
		printf("%g %g %g %.99g %.99g\n", x, y, ulps, my_hypot(x,y), hypot(x, y));
	return 0;
}
#endif