12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485 |
- /* 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 );
- }
|