diff options
Diffstat (limited to 'libm/double/lsqrt.c')
-rw-r--r-- | libm/double/lsqrt.c | 85 |
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 ); +} |