1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
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 );
}
|
