What is the fastest way to find an integer square root using bit shifts?

Question

I was looking for the fastest way to calculate the square root (integer) of a number (integer). I came across this solution on Wikipedia, which finds the square root of a number (if it is a perfect square) or the square root of its nearest lower perfect square (if this number is not a perfect square:

short isqrt(short num) { short res = 0; short bit = 1 << 14; // The second-to-top bit is set: 1L<<30 for long // "bit" starts at the highest power of four <= the argument. while (bit > num) bit >>= 2; while (bit != 0) { if (num >= res + bit) { num -= res + bit; res = (res >> 1) + bit; } else res >>= 1; bit >>= 2; } return res; } 

I tried many test runs to track the algorithm, but it seems I don’t understand the part inside while(bit!=0) . Can anyone explain this part to me?

Answer 1

I also followed a few small examples, and I think I understood. As far as I understand, the algorithm increments the response one binary digit at a time, from the highest bit to the lowest.

Let "num_init" be the num value at the beginning of the function. Suppose that at some iteration we have bit = 4 ^ x, and the number is equal to some value of "num_curr" (a quick look shows that as long as the bit is 0, it is always equal to 4). Then res has the form y * 2 ^ (x + 1), where y ^ 2 + num_curr = num_init, and y is less than the actual answer, but within 2 ^ x.

This invariant with respect to num, res and bit values ​​will be key. The way this is done in code is that

 while (bit != 0) { .... } 

moves our imaginary pointer from left to right, and at each step we determine whether this bit is 0 or 1.

Turning to the first statement of if, suppose that our imaginary "built-in" integer is equal to y, and we look at bit 2 ^ x. Then the bit is 1 if the initial value num is (y + 2 ^ x) ^ 2 = y ^ 2 + y * 2 ^ (x + 1) + 4 ^ x. In other words, a bit is equal to one if the num value at this point is not less than y * 2 ^ (x + 1) + 4 ^ x (since we have the invariant that the num value has fallen by y ^ 2), It is quite convenient res = y * 2 ^ (x + 1) and bit = 4 ^ x. Then we get the point beyond

 if (num >= res + bit) { num -= res + bit; res = (res >> 1) + bit; } else res >>= 1; 

which adds 1 bit in our imaginary place, if necessary, then updates res and num to preserve the invariant. Finally

 bit >>= 2; 

updates the bit and moves everything one step.