summaryrefslogtreecommitdiff
path: root/libm/double/lsqrt.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/lsqrt.c')
-rw-r--r--libm/double/lsqrt.c85
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 );
+}