ReactOS  0.4.15-dev-5455-g015cd25
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  */
33 double
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 
87 int
88 main(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
double _hypot(double x, double y)
Definition: hypot.c:34
GLdouble GLdouble GLdouble r
Definition: gl.h:2055
int main(int argc, char *argv[])
Definition: atactl.cpp:1685
GLdouble GLdouble t
Definition: gl.h:2047
#define DBL_MAX
Definition: gcc_float.h:108
GLint GLint GLint GLint GLint x
Definition: gl.h:1548
#define e
Definition: ke_i.h:82
#define __SQRT_DBL_MAX
Definition: hypot.c:27
#define __SQRT_DBL_MIN
Definition: hypot.c:28
#define printf
Definition: freeldr.h:94
GLdouble GLdouble GLdouble GLdouble q
Definition: gl.h:2063
GLdouble s
Definition: gl.h:2039
_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
GLint GLint GLint GLint GLint GLint y
Definition: gl.h:1548
GLfloat GLfloat p
Definition: glext.h:8902