For floating-point numbers, the square root is almost universally computed by Newton’s method. This method begins with a quadratic function such as f(x)=x2−α to know α since it is a root of f(x). As a short review, Newton’s method computes iteratively to find the approximated root as follows. fn′=xn−xn+1fn⟹xn−xn+1=fn′fn⟹xn+1=xn−fn′fn⟹xn+1=xn−2xnxn2−α=2xnxn2+α=21(xn+xnα)
The iteration converges quadratically, which means that if at some point xn is accurate to n bits, then xn+1 is accurate to 2n bits. So, Newton’s method can be appropriate in terms of integer square root as well. The integer square root here means ⌊α⌋ where α is a 32-bit unsigned integer such that 0≤α≤232−1. Fortunately, Newton’s method works fine in the domain of integers. This is because of the following theorem.
Reasons Why Newton’s Method Works
Let xn+1=⌊21(xn+⌊xnα⌋)⌋. Also, xn and α are positive integers. Then, the following statements hold.
If xn>⌊α⌋, then ⌊α⌋≤xn+1<xn.
If xn=⌊α⌋, then ⌊α⌋≤xn+1≤⌊α⌋+1.
Since xn is an integer, the following holds. xn+1=⌊21(xn+⌊xnα⌋)⌋=⌊21⌊xn+xnα⌋⌋=⌊21(xn+xnα)⌋=⌊2xnxn2+α⌋
The third step is due to the identity discussed in the other post. To prove the first statement, let ϵ be a positive real number such that xn=(1+ϵ)α. Noting that xn>α since xn>⌊α⌋ and xn is an integer, ⌊2xnxn2+α⌋=xn+1≤2xnxn2+α⟹⌊2(1+ϵ)α(1+ϵ)2α+α⌋=xn+1<2xnxn2+xn2⟹⌊2(1+ϵ)ϵ2+2ϵ+2α⌋=xn+1<xn⟹⌊2(1+ϵ)2ϵ+2α⌋≤xn+1<xn⟹⌊α⌋≤xn+1<xn
Now, to prove the second statement, check that α−1<xn≤α since xn=⌊α⌋. So, noting that xn2≤α<(xn+1)2, ⌊2xnxn2+xn2⌋≤⌊2xnxn2+α⌋=xn+1≤⌊2xnxn2+(xn+1)2⌋⟹xn≤xn+1≤⌊xn+1+2xn1⌋⟹⌊α⌋≤xn+1≤xn+1+⌊2xn1⌋⟹⌊α⌋≤xn+1≤xn+1=α+1
As such, these statements imply that if xn is larger than the true value ⌊α⌋, then the next guess xn+1 will be strictly less than the preceding one, but not less than ⌊α⌋. Therefore, the initial guess larger than ⌊α⌋ makes the sequence converge monotonically.
Implementation
The following procedure is the implementation of Newton’s method to find the approximated α. The initial guess x0 is equal to the least power of 2 that is greater than or equal to α. For example, for α=4, x0=2, and for α=5, x0=4.
Because the initial guess is accurate to about 1 bit and Newton’s method converges quadratically, the above code would be expected to converge within about 5 iterations for a 32-bit integer, which requires 4 divisions. It turns out that the maximum number of divisions is 5. If the function to get the number of leading zeros is available, then it is much more simple to set the initial guess.
The following shows the list of initial guesses. Range of x[2,22][22+1,24][24+1,26][26+1,28]⋮[228+1,230][230+1,232−1]Initial Guess2222324⋮215216
Meanwhile, there is a little trick to compute the integer square root for a small integer. In [0,15], the expression (x > 0) + (x > 3) + (x > 8) computes x.
The Integer Square Root of the Integer Square Root
⌊4α⌋=⌊⌊α⌋⌋
It is interesting to note that the integer square root of the integer square root can be also computed from the single integer square root function. To prove the above, let ⌊α⌋=p. Then, p2≤α<(p+1)2 since p is an integer. Also, let ⌊p⌋=q. Then, q2≤p<(q+1)2. So, p+1≤(q+1)2. Therefore, q4≤p2≤α<(p+1)2≤(q+1)4. That is, q4≤α<(q+1)4, which means q≤4α<q+1. Accordingly, ⌊4α⌋=q=⌊p⌋=⌊⌊α⌋⌋.
Reference
[1] Henry S. Warren. 2012. Hacker’s Delight (2nd. ed.). Addison-Wesley Professional.
Keep going!Keep going ×2!Give me more!Thank you, thank youFar too kind!Never gonna give me up?Never gonna let me down?Turn around and desert me!You're an addict!Son of a clapper!No wayGo back to work!This is getting out of handUnbelievablePREPOSTEROUSI N S A N I T YFEED ME A STRAY CAT