1 /* Copyright (C) 1995 DJ Delorie, see COPYING.DJ for details */
3 * hypot() function for DJGPP.
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.
14 * This implementation tries to avoid the above pitfals, without
15 * inflicting too much of a performance hit.
19 #include <msvcrt/float.h>
20 #include <msvcrt/math.h>
21 #include <msvcrt/errno.h>
23 /* Approximate square roots of DBL_MAX and DBL_MIN. Numbers
24 between these two shouldn't neither overflow nor underflow
26 #define __SQRT_DBL_MAX 1.3e+154
27 #define __SQRT_DBL_MIN 2.3e-162
30 _hypot(double x
, double y
)
32 double abig
= fabs(x
), asmall
= fabs(y
);
35 /* Make abig = max(|x|, |y|), asmall = min(|x|, |y|). */
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. */
54 if ((ratio
= asmall
/ abig
) > __SQRT_DBL_MIN
&& abig
< __SQRT_DBL_MAX
)
55 return abig
* sqrt(1.0 + ratio
*ratio
);
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. */
63 double r
= ratio
*ratio
, t
, s
, p
= abig
, q
= asmall
;
72 r
= (q
/ p
) * (q
/ p
);
81 #include <msvcrt/stdio.h>
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
));