home *** CD-ROM | disk | FTP | other *** search
/ PC Extra Super CD 1998 January / PCPLUS131.iso / DJGPP / V2 / DJLSR201.ZIP / src / libc / ansi / math / hypot.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-04  |  3.0 KB  |  101 lines
  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 <math.h>
  21. #include <errno.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
  30. hypot(double x, double y)
  31. {
  32.   double abig = fabs(x), asmall = fabs(y);
  33.   double ratio;
  34.  
  35.   /* Make abig = max(|x|, |y|), asmall = min(|x|, |y|).  */
  36.   if (abig < asmall)
  37.     {
  38.       double temp = abig;
  39.  
  40.       abig = asmall;
  41.       asmall = temp;
  42.     }
  43.  
  44.   /* Trivial case.  */
  45.   if (asmall == 0.)
  46.     return abig;
  47.  
  48.   /* Scale the numbers as much as possible by using its ratio.
  49.      For example, if both ABIG and ASMALL are VERY small, then
  50.      X^2 + Y^2 might be VERY inaccurate due to loss of
  51.      significant digits.  Dividing ASMALL by ABIG scales them
  52.      to a certain degree, so that accuracy is better.  */
  53.  
  54.   if ((ratio = asmall / abig) > __SQRT_DBL_MIN && abig < __SQRT_DBL_MAX)
  55.     return abig * sqrt(1.0 + ratio*ratio);
  56.   else
  57.     {
  58.       /* Slower but safer algorithm due to Moler and Morrison.  Never
  59.          produces any intermediate result greater than roughly the
  60.          larger of X and Y.  Should converge to machine-precision
  61.          accuracy in 3 iterations.  */
  62.  
  63.       double r = ratio*ratio, t, s, p = abig, q = asmall;
  64.  
  65.       do {
  66.         t = 4. + r;
  67.         if (t == 4.)
  68.           break;
  69.         s = r / t;
  70.         p += 2. * s * p;
  71.         q *= s;
  72.         r = (q / p) * (q / p);
  73.       } while (1);
  74.  
  75.       return p;
  76.     }
  77. }
  78.  
  79. #ifdef  TEST
  80.  
  81. #include <stdio.h>
  82.  
  83. int
  84. main(void)
  85. {
  86.   printf("hypot(3, 4) =\t\t\t %25.17e\n", hypot(3., 4.));
  87.   printf("hypot(3*10^150, 4*10^150) =\t %25.17g\n", hypot(3.e+150, 4.e+150));
  88.   printf("hypot(3*10^306, 4*10^306) =\t %25.17g\n", hypot(3.e+306, 4.e+306));
  89.   printf("hypot(3*10^-320, 4*10^-320) =\t %25.17g\n",
  90.          hypot(3.e-320, 4.e-320));
  91.   printf("hypot(0.7*DBL_MAX, 0.7*DBL_MAX) =%25.17g\n",
  92.          hypot(0.7*DBL_MAX, 0.7*DBL_MAX));
  93.   printf("hypot(DBL_MAX, 1.0) =\t\t %25.17g\n", hypot(DBL_MAX, 1.0));
  94.   printf("hypot(1.0, DBL_MAX) =\t\t %25.17g\n", hypot(1.0, DBL_MAX));
  95.   printf("hypot(0.0, DBL_MAX) =\t\t %25.17g\n", hypot(0.0, DBL_MAX));
  96.  
  97.   return 0;
  98. }
  99.  
  100. #endif
  101.