lib/crtdll/math/hypot.c

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