ReactOS 0.4.15-dev-8219-ge8b88cf
hypot.c
Go to the documentation of this file.
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#include <precomp.h>
19
20#if (_MSC_VER >= 1920)
21#pragma function(_hypot)
22#endif
23
24/* Approximate square roots of DBL_MAX and DBL_MIN. Numbers
25 between these two shouldn't neither overflow nor underflow
26 when squared. */
27#define __SQRT_DBL_MAX 1.3e+154
28#define __SQRT_DBL_MIN 2.3e-162
29
30/*
31 * @implemented
32 */
33double
34_hypot(double x, double y)
35{
36 double abig = fabs(x), asmall = fabs(y);
37 double ratio;
38
39 /* Make abig = max(|x|, |y|), asmall = min(|x|, |y|). */
40 if (abig < asmall)
41 {
42 double temp = abig;
43
44 abig = asmall;
45 asmall = temp;
46 }
47
48 /* Trivial case. */
49 if (asmall == 0.)
50 return abig;
51
52 /* Scale the numbers as much as possible by using its ratio.
53 For example, if both ABIG and ASMALL are VERY small, then
54 X^2 + Y^2 might be VERY inaccurate due to loss of
55 significant digits. Dividing ASMALL by ABIG scales them
56 to a certain degree, so that accuracy is better. */
57
58 if ((ratio = asmall / abig) > __SQRT_DBL_MIN && abig < __SQRT_DBL_MAX)
59 return abig * sqrt(1.0 + ratio*ratio);
60 else
61 {
62 /* Slower but safer algorithm due to Moler and Morrison. Never
63 produces any intermediate result greater than roughly the
64 larger of X and Y. Should converge to machine-precision
65 accuracy in 3 iterations. */
66
67 double r = ratio*ratio, t, s, p = abig, q = asmall;
68
69 do {
70 t = 4. + r;
71 if (t == 4.)
72 break;
73 s = r / t;
74 p += 2. * s * p;
75 q *= s;
76 r = (q / p) * (q / p);
77 } while (1);
78
79 return p;
80 }
81}
82
83#ifdef TEST
84
85#include <msvcrt/stdio.h>
86
87int
88main(void)
89{
90 printf("hypot(3, 4) =\t\t\t %25.17e\n", _hypot(3., 4.));
91 printf("hypot(3*10^150, 4*10^150) =\t %25.17g\n", _hypot(3.e+150, 4.e+150));
92 printf("hypot(3*10^306, 4*10^306) =\t %25.17g\n", _hypot(3.e+306, 4.e+306));
93 printf("hypot(3*10^-320, 4*10^-320) =\t %25.17g\n",_hypot(3.e-320, 4.e-320));
94 printf("hypot(0.7*DBL_MAX, 0.7*DBL_MAX) =%25.17g\n",_hypot(0.7*DBL_MAX, 0.7*DBL_MAX));
95 printf("hypot(DBL_MAX, 1.0) =\t\t %25.17g\n", _hypot(DBL_MAX, 1.0));
96 printf("hypot(1.0, DBL_MAX) =\t\t %25.17g\n", _hypot(1.0, DBL_MAX));
97 printf("hypot(0.0, DBL_MAX) =\t\t %25.17g\n", _hypot(0.0, DBL_MAX));
98
99 return 0;
100}
101
102#endif
_STLP_DECLSPEC complex< float > _STLP_CALL sqrt(const complex< float > &)
Definition: complex.cpp:188
int main()
Definition: test.c:6
#define printf
Definition: freeldr.h:97
#define DBL_MAX
Definition: gcc_float.h:108
GLint GLint GLint GLint GLint x
Definition: gl.h:1548
GLdouble s
Definition: gl.h:2039
GLint GLint GLint GLint GLint GLint y
Definition: gl.h:1548
GLdouble GLdouble GLdouble r
Definition: gl.h:2055
GLdouble GLdouble t
Definition: gl.h:2047
GLdouble GLdouble GLdouble GLdouble q
Definition: gl.h:2063
GLfloat GLfloat p
Definition: glext.h:8902
#define __SQRT_DBL_MAX
Definition: hypot.c:27
double _hypot(double x, double y)
Definition: hypot.c:34
#define __SQRT_DBL_MIN
Definition: hypot.c:28
_Check_return_ _CRT_JIT_INTRINSIC double __cdecl fabs(_In_ double x)
Definition: fabs.c:17
static calc_node_t temp
Definition: rpn_ieee.c:38