/* hypot - sqrt(x * x + y * y)
**
** Because this was omitted from the ANSI standards, this version
** is included for those systems that do not include hypot as an
** extension to libm.a. Note: GNU version was not used because it
** was not properly coded to minimize potential overflow.
**
** The proper technique for determining hypot is to factor out the
** larger of the two terms, thus leaving a possible case of float
** overflow when max(x,y)*sqrt(2) > max machine value. This allows
** a wider range of numbers than the alternative of the sum of the
** squares < max machine value. For an Intel x87 IEEE double of
** approximately 1.8e308, only argument values > 1.27e308 are at
** risk of causing overflow. Whereas, not using this method limits
** the range to values less that 9.5e153 --- a considerable reduction
** in range!
*/
extern double sqrt(double);
double
hypot(double x, double y) {
if ( x < 0.)
x = -x;
else if (x == 0.)
return (y < 0. ? -y : y);
if (y < 0.)
y = -y;
else if (y == 0.)
return (x);
if ( x < y ) {
x /= y;
return ( y * sqrt( 1. + x * x ) );
} else {
y /= x;
return ( x * sqrt( 1. + y * y ) );
}
}