1 /* Copyright (C) 1995 DJ Delorie, see COPYING.DJ for details */ 2 /* 3 * hypot() function for DJGPP. 4 * 5 * hypot() computes sqrt(x^2 + y^2). The problem with the obvious 6 * naive implementation is that it might fail for very large or 7 * very small arguments. For instance, for large x or y the result 8 * might overflow even if the value of the function should not, 9 * because squaring a large number might trigger an overflow. For 10 * very small numbers, their square might underflow and will be 11 * silently replaced by zero; this won't cause an exception, but might 12 * have an adverse effect on the accuracy of the result. 13 * 14 * This implementation tries to avoid the above pitfals, without 15 * inflicting too much of a performance hit. 16 * 17 */ 18 #include <precomp.h> 19 20 #if (_MSC_VER >= 1920) 21 #pragma function(_hypot) 22 #endif 23 24 /* Approximate square roots of DBL_MAX and DBL_MIN. Numbers 25 between these two shouldn't neither overflow nor underflow 26 when squared. */ 27 #define __SQRT_DBL_MAX 1.3e+154 28 #define __SQRT_DBL_MIN 2.3e-162 29 30 /* 31 * @implemented 32 */ 33 double 34 _hypot(double x, double y) 35 { 36 double abig = fabs(x), asmall = fabs(y); 37 double ratio; 38 39 /* Make abig = max(|x|, |y|), asmall = min(|x|, |y|). */ 40 if (abig < asmall) 41 { 42 double temp = abig; 43 44 abig = asmall; 45 asmall = temp; 46 } 47 48 /* Trivial case. */ 49 if (asmall == 0.) 50 return abig; 51 52 /* Scale the numbers as much as possible by using its ratio. 53 For example, if both ABIG and ASMALL are VERY small, then 54 X^2 + Y^2 might be VERY inaccurate due to loss of 55 significant digits. Dividing ASMALL by ABIG scales them 56 to a certain degree, so that accuracy is better. */ 57 58 if ((ratio = asmall / abig) > __SQRT_DBL_MIN && abig < __SQRT_DBL_MAX) 59 return abig * sqrt(1.0 + ratio*ratio); 60 else 61 { 62 /* Slower but safer algorithm due to Moler and Morrison. Never 63 produces any intermediate result greater than roughly the 64 larger of X and Y. Should converge to machine-precision 65 accuracy in 3 iterations. */ 66 67 double r = ratio*ratio, t, s, p = abig, q = asmall; 68 69 do { 70 t = 4. + r; 71 if (t == 4.) 72 break; 73 s = r / t; 74 p += 2. * s * p; 75 q *= s; 76 r = (q / p) * (q / p); 77 } while (1); 78 79 return p; 80 } 81 } 82 83 #ifdef TEST 84 85 #include <msvcrt/stdio.h> 86 87 int 88 main(void) 89 { 90 printf("hypot(3, 4) =\t\t\t %25.17e\n", _hypot(3., 4.)); 91 printf("hypot(3*10^150, 4*10^150) =\t %25.17g\n", _hypot(3.e+150, 4.e+150)); 92 printf("hypot(3*10^306, 4*10^306) =\t %25.17g\n", _hypot(3.e+306, 4.e+306)); 93 printf("hypot(3*10^-320, 4*10^-320) =\t %25.17g\n",_hypot(3.e-320, 4.e-320)); 94 printf("hypot(0.7*DBL_MAX, 0.7*DBL_MAX) =%25.17g\n",_hypot(0.7*DBL_MAX, 0.7*DBL_MAX)); 95 printf("hypot(DBL_MAX, 1.0) =\t\t %25.17g\n", _hypot(DBL_MAX, 1.0)); 96 printf("hypot(1.0, DBL_MAX) =\t\t %25.17g\n", _hypot(1.0, DBL_MAX)); 97 printf("hypot(0.0, DBL_MAX) =\t\t %25.17g\n", _hypot(0.0, DBL_MAX)); 98 99 return 0; 100 } 101 102 #endif 103