# Subversion RepositoriesKolibri OS

Rev
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.