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