forked from KolibriOS/kolibrios
101 lines
2.9 KiB
C
101 lines
2.9 KiB
C
|
/* Copyright (C) 1995 DJ Delorie, see COPYING.DJ for details */
|
||
|
|
/*
|
||
|
|
* hypot() function for DJGPP.
|
||
|
|
*
|
||
|
|
* hypot() computes sqrt(x^2 + y^2). The problem with the obvious
|
||
|
|
* naive implementation is that it might fail for very large or
|
||
|
|
* very small arguments. For instance, for large x or y the result
|
||
|
|
* might overflow even if the value of the function should not,
|
||
|
|
* because squaring a large number might trigger an overflow. For
|
||
|
|
* very small numbers, their square might underflow and will be
|
||
|
|
* silently replaced by zero; this won't cause an exception, but might
|
||
|
|
* have an adverse effect on the accuracy of the result.
|
||
|
|
*
|
||
|
|
* This implementation tries to avoid the above pitfals, without
|
||
|
|
* inflicting too much of a performance hit.
|
||
|
|
*
|
||
|
|
*/
|
||
|
|
|
||
|
|
/// #include <float.h>
|
||
|
|
#include <math.h>
|
||
|
|
#include <errno.h>
|
||
|
|
|
||
|
|
/* Approximate square roots of DBL_MAX and DBL_MIN. Numbers
|
||
|
|
between these two shouldn't neither overflow nor underflow
|
||
|
|
when squared. */
|
||
|
|
#define __SQRT_DBL_MAX 1.3e+154
|
||
|
|
#define __SQRT_DBL_MIN 2.3e-162
|
||
|
|
|
||
|
|
double
|
||
|
|
hypot(double x, double y)
|
||
|
|
{
|
||
|
|
double abig = fabs(x), asmall = fabs(y);
|
||
|
|
double ratio;
|
||
|
|
|
||
|
|
/* Make abig = max(|x|, |y|), asmall = min(|x|, |y|). */
|
||
|
|
if (abig < asmall)
|
||
|
|
{
|
||
|
|
double temp = abig;
|
||
|
|
|
||
|
|
abig = asmall;
|
||
|
|
asmall = temp;
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Trivial case. */
|
||
|
|
if (asmall == 0.)
|
||
|
|
return abig;
|
||
|
|
|
||
|
|
/* Scale the numbers as much as possible by using its ratio.
|
||
|
|
For example, if both ABIG and ASMALL are VERY small, then
|
||
|
|
X^2 + Y^2 might be VERY inaccurate due to loss of
|
||
|
|
significant digits. Dividing ASMALL by ABIG scales them
|
||
|
|
to a certain degree, so that accuracy is better. */
|
||
|
|
|
||
|
|
if ((ratio = asmall / abig) > __SQRT_DBL_MIN && abig < __SQRT_DBL_MAX)
|
||
|
|
return abig * sqrt(1.0 + ratio*ratio);
|
||
|
|
else
|
||
|
|
{
|
||
|
|
/* Slower but safer algorithm due to Moler and Morrison. Never
|
||
|
|
produces any intermediate result greater than roughly the
|
||
|
|
larger of X and Y. Should converge to machine-precision
|
||
|
|
accuracy in 3 iterations. */
|
||
|
|
|
||
|
|
double r = ratio*ratio, t, s, p = abig, q = asmall;
|
||
|
|
|
||
|
|
do {
|
||
|
|
t = 4. + r;
|
||
|
|
if (t == 4.)
|
||
|
|
break;
|
||
|
|
s = r / t;
|
||
|
|
p += 2. * s * p;
|
||
|
|
q *= s;
|
||
|
|
r = (q / p) * (q / p);
|
||
|
|
} while (1);
|
||
|
|
|
||
|
|
return p;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
#ifdef TEST
|
||
|
|
|
||
|
|
#include <stdio.h>
|
||
|
|
|
||
|
|
int
|
||
|
|
main(void)
|
||
|
|
{
|
||
|
|
printf("hypot(3, 4) =\t\t\t %25.17e\n", hypot(3., 4.));
|
||
|
|
printf("hypot(3*10^150, 4*10^150) =\t %25.17g\n", hypot(3.e+150, 4.e+150));
|
||
|
|
printf("hypot(3*10^306, 4*10^306) =\t %25.17g\n", hypot(3.e+306, 4.e+306));
|
||
|
|
printf("hypot(3*10^-320, 4*10^-320) =\t %25.17g\n",
|
||
|
|
hypot(3.e-320, 4.e-320));
|
||
|
|
printf("hypot(0.7*DBL_MAX, 0.7*DBL_MAX) =%25.17g\n",
|
||
|
|
hypot(0.7*DBL_MAX, 0.7*DBL_MAX));
|
||
|
|
printf("hypot(DBL_MAX, 1.0) =\t\t %25.17g\n", hypot(DBL_MAX, 1.0));
|
||
|
|
printf("hypot(1.0, DBL_MAX) =\t\t %25.17g\n", hypot(1.0, DBL_MAX));
|
||
|
|
printf("hypot(0.0, DBL_MAX) =\t\t %25.17g\n", hypot(0.0, DBL_MAX));
|
||
|
|
|
||
|
|
return 0;
|
||
|
|
}
|
||
|
|
|
||
|
|
#endif
|