ReactOS  0.4.13-dev-242-g611e6d7
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 /* Approximate square roots of DBL_MAX and DBL_MIN. Numbers
21  between these two shouldn't neither overflow nor underflow
22  when squared. */
23 #define __SQRT_DBL_MAX 1.3e+154
24 #define __SQRT_DBL_MIN 2.3e-162
25 
26 /*
27  * @implemented
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 <msvcrt/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
_STLP_DECLSPEC complex< float > _STLP_CALL sqrt(const complex< float > &)
Definition: complex.cpp:188
double _hypot(double x, double y)
Definition: hypot.c:30
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:23
#define __SQRT_DBL_MIN
Definition: hypot.c:24
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)
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
#define printf
Definition: config.h:203