uclibc-ng/libm/double/lsqrt.c

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