1 files changed, 85 insertions, 0 deletions
diff --git a/libm/double/lsqrt.c b/libm/double/lsqrt.c
new file mode 100644
index 000000000..bf85a54f1
--- /dev/null
+++ b/libm/double/lsqrt.c
@@ -0,0 +1,85 @@
+/*							lsqrt.c
+ *
+ *	Integer square root
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long x, y;
+ * long lsqrt();
+ *
+ * y = lsqrt( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns a long integer square root of the long integer
+ * argument.  The computation is by binary long division.
+ *
+ * The largest possible result is lsqrt(2,147,483,647)
+ * = 46341.
+ *
+ * If x < 0, the square root of |x| is returned, and an
+ * error message is printed.
+ *
+ *
+ * ACCURACY:
+ *
+ * An extra, roundoff, bit is computed; hence the result
+ * is the nearest integer to the actual square root.
+ * NOTE: only DEC arithmetic is currently supported.
+ *
+ */
+
+/*
+Cephes Math Library Release 2.0:  April, 1987
+Copyright 1984, 1987 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+#include <math.h>
+
+long lsqrt(x)
+long x;
+{
+long num, sq;
+long temp;
+int i, j, k, n;
+
+if( x < 0 )
+	{
+	mtherr( "lsqrt", DOMAIN );
+	x = -x;
+	}
+
+num = 0;
+sq = 0;
+k = 24;
+n = 4;
+
+for( j=0; j<4; j++ )
+	{
+	num |= (x >> k) & 0xff;	/* bring in next byte of arg */
+	if( j == 3 )		/* do roundoff bit at end */
+		n = 5;
+	for( i=0; i<n; i++ )
+		{
+		num <<= 2;		/* next 2 bits of arg */
+		sq <<= 1;		/* shift up answer */
+		temp = (sq << 1) + 256;	/* trial divisor */
+		temp = num - temp;
+		if( temp >= 0 )
+			{
+			num = temp;	/* it went in */
+			sq += 256;	/* answer bit = 1 */
+			}
+		}
+	k -= 8;	/* shift count to get next byte of arg */
+	}
+
+sq += 256;	/* add roundoff bit */
+sq >>= 9;	/* truncate */
+return( sq );
+}