Subversion Repositories Kolibri OS

Rev

Rev 8793 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed

  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.  
  19. /// #include <float.h>
  20. #include <errno.h>
  21. #include <math.h>
  22.  
  23. /* Approximate square roots of DBL_MAX and DBL_MIN.  Numbers
  24.    between these two shouldn't neither overflow nor underflow
  25.    when squared.  */
  26. #define __SQRT_DBL_MAX 1.3e+154
  27. #define __SQRT_DBL_MIN 2.3e-162
  28.  
  29. double hypot(double x, double y)
  30. {
  31.     double abig = fabs(x), asmall = fabs(y);
  32.     double ratio;
  33.  
  34.     /* Make abig = max(|x|, |y|), asmall = min(|x|, |y|).  */
  35.     if (abig < asmall) {
  36.         double temp = abig;
  37.  
  38.         abig = asmall;
  39.         asmall = temp;
  40.     }
  41.  
  42.     /* Trivial case.  */
  43.     if (asmall == 0.)
  44.         return abig;
  45.  
  46.     /* Scale the numbers as much as possible by using its ratio.
  47.        For example, if both ABIG and ASMALL are VERY small, then
  48.        X^2 + Y^2 might be VERY inaccurate due to loss of
  49.        significant digits.  Dividing ASMALL by ABIG scales them
  50.        to a certain degree, so that accuracy is better.  */
  51.  
  52.     if ((ratio = asmall / abig) > __SQRT_DBL_MIN && abig < __SQRT_DBL_MAX)
  53.         return abig * sqrt(1.0 + ratio * ratio);
  54.     else {
  55.         /* Slower but safer algorithm due to Moler and Morrison.  Never
  56.            produces any intermediate result greater than roughly the
  57.            larger of X and Y.  Should converge to machine-precision
  58.            accuracy in 3 iterations.  */
  59.  
  60.         double r = ratio * ratio, t, s, p = abig, q = asmall;
  61.  
  62.         do {
  63.             t = 4. + r;
  64.             if (t == 4.)
  65.                 break;
  66.             s = r / t;
  67.             p += 2. * s * p;
  68.             q *= s;
  69.             r = (q / p) * (q / p);
  70.         } while (1);
  71.  
  72.         return p;
  73.     }
  74. }
  75.  
  76. #ifdef TEST
  77.  
  78. #include <stdio.h>
  79.  
  80. int main(void)
  81. {
  82.     printf("hypot(3, 4) =\t\t\t %25.17e\n", hypot(3., 4.));
  83.     printf("hypot(3*10^150, 4*10^150) =\t %25.17g\n", hypot(3.e+150, 4.e+150));
  84.     printf("hypot(3*10^306, 4*10^306) =\t %25.17g\n", hypot(3.e+306, 4.e+306));
  85.     printf("hypot(3*10^-320, 4*10^-320) =\t %25.17g\n",
  86.         hypot(3.e-320, 4.e-320));
  87.     printf("hypot(0.7*DBL_MAX, 0.7*DBL_MAX) =%25.17g\n",
  88.         hypot(0.7 * DBL_MAX, 0.7 * DBL_MAX));
  89.     printf("hypot(DBL_MAX, 1.0) =\t\t %25.17g\n", hypot(DBL_MAX, 1.0));
  90.     printf("hypot(1.0, DBL_MAX) =\t\t %25.17g\n", hypot(1.0, DBL_MAX));
  91.     printf("hypot(0.0, DBL_MAX) =\t\t %25.17g\n", hypot(0.0, DBL_MAX));
  92.  
  93.     return 0;
  94. }
  95.  
  96. #endif
  97.