lib/msvcrt/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 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 <crtdll/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",_hypot(3.e-320, 4.e-320));
  90   printf("hypot(0.7*DBL_MAX, 0.7*DBL_MAX) =%25.17g\n",_hypot(0.7*DBL_MAX, 0.7*DBL_MAX));
  91   printf("hypot(DBL_MAX, 1.0) =\t\t %25.17g\n", _hypot(DBL_MAX, 1.0));
  92   printf("hypot(1.0, DBL_MAX) =\t\t %25.17g\n", _hypot(1.0, DBL_MAX));
  93   printf("hypot(0.0, DBL_MAX) =\t\t %25.17g\n", _hypot(0.0, DBL_MAX));
  94
  95   return 0;
  96 }
  97
  98 #endif