/* Copyright (C) 1995 DJ Delorie, see COPYING.DJ for details */ /* * hypot() function for DJGPP. * * hypot() computes sqrt(x^2 + y^2). The problem with the obvious * naive implementation is that it might fail for very large or * very small arguments. For instance, for large x or y the result * might overflow even if the value of the function should not, * because squaring a large number might trigger an overflow. For * very small numbers, their square might underflow and will be * silently replaced by zero; this won't cause an exception, but might * have an adverse effect on the accuracy of the result. * * This implementation tries to avoid the above pitfals, without * inflicting too much of a performance hit. * */ /// #include <float.h> #include <math.h> #include <errno.h> /* Approximate square roots of DBL_MAX and DBL_MIN. Numbers between these two shouldn't neither overflow nor underflow when squared. */ #define __SQRT_DBL_MAX 1.3e+154 #define __SQRT_DBL_MIN 2.3e-162 double hypot(double x, double y) { double abig = fabs(x), asmall = fabs(y); double ratio; /* Make abig = max(|x|, |y|), asmall = min(|x|, |y|). */ if (abig < asmall) { double temp = abig; abig = asmall; asmall = temp; } /* Trivial case. */ if (asmall == 0.) return abig; /* Scale the numbers as much as possible by using its ratio. For example, if both ABIG and ASMALL are VERY small, then X^2 + Y^2 might be VERY inaccurate due to loss of significant digits. Dividing ASMALL by ABIG scales them to a certain degree, so that accuracy is better. */ if ((ratio = asmall / abig) > __SQRT_DBL_MIN && abig < __SQRT_DBL_MAX) return abig * sqrt(1.0 + ratio*ratio); else { /* Slower but safer algorithm due to Moler and Morrison. Never produces any intermediate result greater than roughly the larger of X and Y. Should converge to machine-precision accuracy in 3 iterations. */ double r = ratio*ratio, t, s, p = abig, q = asmall; do { t = 4. + r; if (t == 4.) break; s = r / t; p += 2. * s * p; q *= s; r = (q / p) * (q / p); } while (1); return p; } } #ifdef TEST #include <stdio.h> int main(void) { printf("hypot(3, 4) =\t\t\t %25.17e\n", hypot(3., 4.)); printf("hypot(3*10^150, 4*10^150) =\t %25.17g\n", hypot(3.e+150, 4.e+150)); printf("hypot(3*10^306, 4*10^306) =\t %25.17g\n", hypot(3.e+306, 4.e+306)); printf("hypot(3*10^-320, 4*10^-320) =\t %25.17g\n", hypot(3.e-320, 4.e-320)); printf("hypot(0.7*DBL_MAX, 0.7*DBL_MAX) =%25.17g\n", hypot(0.7*DBL_MAX, 0.7*DBL_MAX)); printf("hypot(DBL_MAX, 1.0) =\t\t %25.17g\n", hypot(DBL_MAX, 1.0)); printf("hypot(1.0, DBL_MAX) =\t\t %25.17g\n", hypot(1.0, DBL_MAX)); printf("hypot(0.0, DBL_MAX) =\t\t %25.17g\n", hypot(0.0, DBL_MAX)); return 0; } #endif