1 files changed, 39 insertions, 0 deletions
diff --git a/src/hypot.c b/src/hypot.c
new file mode 100644
index 00000000..269c146b
--- /dev/null
+++ b/src/hypot.c
@@ -0,0 +1,39 @@
+#ifndef lint
+static const char SCCSID[]="@(#)hypot.c	4.4	93/06/12	GIE	REL";
+#endif
+/* 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 ) );
+	}
+}