Rev 193 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
259 | daniel-mar | 1 | // From: Ron_Hunsinger@bmug.org (Ron Hunsinger) |
2 | |||
3 | #define BITS 32 |
||
4 | // BITS is the smallest even integer such that 2^BITS > ULONG_MAX. |
||
5 | // If you port this routine to a computer where 2^(BITS-1) > ULONG_MAX |
||
6 | // (Unisys A-series, for example, where integers would generally be |
||
7 | // stored in the 39-bit mantissa of a 48-bit word), unwind the first |
||
8 | // iteration of the loop, and special-case references to 'half' in that |
||
9 | // first iteration so as to avoid overflow. |
||
10 | |||
11 | unsigned long lsqrt (unsigned long n) { |
||
12 | // Compute the integer square root of the integer argument n |
||
13 | // Method is to divide n by x computing the quotient x and remainder r |
||
14 | // Notice that the divisor x is changing as the quotient x changes |
||
15 | |||
16 | // Instead of shifting the dividend/remainder left, we shift the |
||
17 | // quotient/divisor right. The binary point starts at the extreme |
||
18 | // left, and shifts two bits at a time to the extreme right. |
||
19 | |||
20 | // The residue contains n-x^2. (Within these comments, the ^ operator |
||
21 | // signifies exponentiation rather than exclusive or. Also, the / |
||
22 | // operator returns fractions, rather than truncating, so 1/4 means |
||
23 | // one fourth, not zero.) |
||
24 | |||
25 | // Since (x + 1/2)^2 == x^2 + x + 1/4, |
||
26 | // n - (x + 1/2)^2 == (n - x^2) - (x + 1/4) |
||
27 | // Thus, we can increase x by 1/2 if we decrease (n-x^2) by (x+1/4) |
||
28 | |||
29 | unsigned long residue; // n - x^2 |
||
30 | unsigned long root; // x + 1/4 |
||
31 | unsigned long half; // 1/2 |
||
32 | |||
33 | residue = n; // n - (x = 0)^2, with suitable alignment |
||
34 | root = 1 << (BITS - 2); // x + 1/4, shifted left BITS places |
||
35 | half = root + root; // 1/2, shifted left BITS places |
||
36 | |||
37 | do { |
||
38 | if (root <= residue) { // Whenever we can, |
||
39 | residue -= root; // decrease (n-x^2) by (x+1/4) |
||
40 | root += half; } // increase x by 1/2 |
||
41 | half >>= 2; // Shift binary point 2 places right |
||
42 | root -= half; // x{+1/2}+1/4 - 1/8 == x{+1/2}+1/8 |
||
43 | root >>= 1; // 2x{+1}+1/4, shifted right 2 places |
||
44 | } while (half); // When 1/2 == 0, binary point is at far right |
||
45 | |||
46 | // Omit the following line to truncate instead of rounding |
||
47 | if (root < residue) ++root; |
||
48 | |||
49 | return root; // Guaranteed to be correctly rounded (or truncated) |
||
50 | } |