Rev 193 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 193 | Rev 259 | ||
---|---|---|---|
Line 1... | Line -... | ||
1 | // From: Ron_Hunsinger@bmug.org (Ron Hunsinger) |
- | |
2 | 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 | } |
|
- | 51 |