Subversion Repositories filter_foundry

#define BITS 32

// BITS is the smallest even integer such that 2^BITS > ULONG_MAX.

// If you port this routine to a computer where 2^(BITS-1) > ULONG_MAX

// (Unisys A-series, for example, where integers would generally be

// stored in the 39-bit mantissa of a 48-bit word), unwind the first 

// iteration of the loop, and special-case references to 'half' in that 

// first iteration so as to avoid overflow.

unsigned long lsqrt (unsigned long n) {

    // Compute the integer square root of the integer argument n

    // Method is to divide n by x computing the quotient x and remainder r

    // Notice that the divisor x is changing as the quotient x changes

    // Instead of shifting the dividend/remainder left, we shift the

    // quotient/divisor right.  The binary point starts at the extreme

    // left, and shifts two bits at a time to the extreme right.

    // The residue contains n-x^2.  (Within these comments, the ^ operator

    // signifies exponentiation rather than exclusive or.  Also, the /

    // operator returns fractions, rather than truncating, so 1/4 means

    // one fourth, not zero.)

    // Since (x + 1/2)^2 == x^2 + x + 1/4,

    //   n - (x + 1/2)^2 == (n - x^2) - (x + 1/4)

    // Thus, we can increase x by 1/2 if we decrease (n-x^2) by (x+1/4)

    unsigned long residue;   // n - x^2

    unsigned long root;      // x + 1/4

    unsigned long half;      // 1/2

    residue = n;             // n - (x = 0)^2, with suitable alignment

    root = 1 << (BITS - 2);  // x + 1/4, shifted left BITS places

    half = root + root;      // 1/2, shifted left BITS places

    do {

        if (root <= residue) {   // Whenever we can,

            residue -= root;         // decrease (n-x^2) by (x+1/4)

            root += half; }          // increase x by 1/2

        half >>= 2;          // Shift binary point 2 places right

        root -= half;        // x{+1/2}+1/4 - 1/8 == x{+1/2}+1/8

        root >>= 1;          // 2x{+1}+1/4, shifted right 2 places

        } while (half);      // When 1/2 == 0, binary point is at far right

    // Omit the following line to truncate instead of rounding

    if (root < residue) ++root;

    return root;   // Guaranteed to be correctly rounded (or truncated)

    }