Integer Square Root

Newton’s Method

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) = x^2 - \alpha$ to know $\sqrt{\alpha}$ since it is a root of $f(x)$. As a short review, Newton’s method computes iteratively to find the approximated root as follows. \begin{aligned} f'_n = \frac{f_n}{x_n - x_{n + 1}} &\implies x_n - x_{n + 1} = \frac{f_n}{f'_n} \implies x_{n + 1} = x_n - \frac{f_n}{f'_n} \\\\ &\implies x_{n + 1} = x_n - \frac{x_n^2 - \alpha}{2 x_n} = \frac{x_n^2 + \alpha}{2 x_n} = \frac{1}{2} \left( x_n + \frac{\alpha}{x_n} \right) \end{aligned}

The iteration converges quadratically, which means that if at some point $x_n$ is accurate to $n$ bits, then $x_{n + 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 $\lfloor \sqrt{\alpha} \rfloor$ where $\alpha$ is a 32-bit unsigned integer such that $0 \leq \alpha \leq 2^{32} - 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 $x_{n + 1} = \lfloor \frac{1}{2} ( x_n + \lfloor \frac{\alpha}{x_n} \rfloor ) \rfloor$. Also, $x_n$ and $\alpha$ are positive integers. Then, the following statements hold.

• If $x_n > \lfloor \sqrt{\alpha} \rfloor$, then $\lfloor \sqrt{\alpha} \rfloor \leq x_{n + 1} < x_n$.
• If $x_n = \lfloor \sqrt{\alpha} \rfloor$, then $\lfloor \sqrt{\alpha} \rfloor \leq x_{n + 1} \leq \lfloor \sqrt{\alpha} \rfloor + 1$.

Since $x_n$ is an integer, the following holds. \begin{aligned} x_{n + 1} = \left \lfloor \cfrac{1}{2} \left( x_n + \left \lfloor \cfrac{\alpha}{x_n} \right \rfloor \right) \right \rfloor = \left \lfloor \cfrac{1}{2} \left \lfloor x_n + \cfrac{\alpha}{x_n} \right \rfloor \right \rfloor = \left \lfloor \cfrac{1}{2} \left( x_n + \cfrac{\alpha}{x_n} \right) \right \rfloor = \left \lfloor \cfrac{x_n^2 + \alpha}{2 x_n} \right \rfloor \end{aligned}

The third step is due to the identity discussed in the other post. To prove the first statement, let $\epsilon$ be a positive real number such that $x_n = (1 + \epsilon) \sqrt{\alpha}$. Noting that $x_n > \sqrt{\alpha}$ since $x_n > \lfloor \sqrt{\alpha} \rfloor$ and $x_n$ is an integer, \begin{aligned} \left \lfloor \cfrac{x_n^2 + \alpha}{2 x_n} \right \rfloor = x_{n + 1} \leq \cfrac{x_n^2 + \alpha}{2 x_n} &\implies \left \lfloor \cfrac{(1 + \epsilon)^2 \alpha + \alpha}{2 (1 + \epsilon) \sqrt{\alpha}} \right \rfloor = x_{n + 1} < \cfrac{x_n^2 + x_n^2}{2 x_n} \\ &\implies \left \lfloor \cfrac{\epsilon^2 + 2 \epsilon + 2}{2 (1 + \epsilon)} \sqrt{\alpha} \right \rfloor = x_{n + 1} < x_n \\ &\implies \left \lfloor \cfrac{2 \epsilon + 2}{2 (1 + \epsilon)} \sqrt{\alpha} \right \rfloor \leq x_{n + 1} < x_n \\ &\implies \left \lfloor \sqrt{\alpha} \right \rfloor \leq x_{n + 1} < x_n \end{aligned}

Now, to prove the second statement, check that $\sqrt{\alpha} - 1 < x_n \leq \sqrt{\alpha}$ since $x_n = \lfloor \sqrt{\alpha} \rfloor$. So, noting that $x_n^2 \leq \alpha < (x_n + 1)^2$, \begin{aligned} \left \lfloor \cfrac{x_n^2 + x_n^2}{2 x_n} \right \rfloor \leq \left \lfloor \cfrac{x_n^2 + \alpha}{2 x_n} \right \rfloor = x_{n + 1} \leq \left \lfloor \cfrac{x_n^2 + (x_n + 1)^2}{2 x_n} \right \rfloor &\implies x_n \leq x_{n + 1} \leq \left \lfloor x_n + 1 + \cfrac{1}{2 x_n} \right \rfloor \\ &\implies \left \lfloor \sqrt{\alpha} \right \rfloor \leq x_{n + 1} \leq x_n + 1 + \left \lfloor \cfrac{1}{2 x_n} \right \rfloor\\ &\implies \left \lfloor \sqrt{\alpha} \right \rfloor \leq x_{n + 1} \leq x_n + 1 = \sqrt{\alpha} + 1 \end{aligned}

As such, these statements imply that if $x_n$ is larger than the true value $\lfloor \sqrt{\alpha} \rfloor$, then the next guess $x_{n + 1}$ will be strictly less than the preceding one, but not less than $\lfloor \sqrt{\alpha} \rfloor$. Therefore, the initial guess larger than $\lfloor \sqrt{\alpha} \rfloor$ makes the sequence converge monotonically.

Implementation

The following procedure is the implementation of Newton’s method to find the approximated $\sqrt{\alpha}$. The initial guess $x_0$ is equal to the least power of 2 that is greater than or equal to $\sqrt{\alpha}$. For example, for $\alpha = 4$, $x_0 = 2$, and for $\alpha = 5$, $x_0 = 4$.

int sqrtInteger(unsigned int a)
{
if (a <= 1) return a;

int n = 1;
unsigned int guess = a - 1;
if (guess > 65535) { n += 8; guess >>= 16; }
if (guess > 255) { n += 4; guess >>= 8; }
if (guess > 15) { n += 2; guess >>= 4; }
if (guess > 3) n++;

unsigned int x_curr = 1 << n;
unsigned int x_next = (x_curr + (a >> n)) >> 1;
while (x_next < x_curr) {
x_curr = x_next;
x_next = (x_curr + (a / x_curr)) >> 1;
}
return x_curr;
}


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.

int sqrtInteger(unsigned int a)
{
if (a <= 1) return a;

int n = 16 - clz( a - 1 ) / 2;
unsigned int x_curr = 1 << n;
unsigned int x_next = (x_curr + (a >> n)) >> 1;
while (x_next < x_curr) {
x_curr = x_next;
x_next = (x_curr + (a / x_curr)) >> 1;
}
return x_curr;
}


The following shows the list of initial guesses. \begin{aligned} \begin{array}{ccc} \text{Range of } x & \text{Initial Guess} \\ \hline \\ [2, 2^2] & 2 \\\\ [2^2 + 1, 2^4] & 2^2 \\\\ [2^4 + 1, 2^6] & 2^3 \\\\ [2^6 + 1, 2^8] & 2^4 \\\\ \vdots & \vdots \\\\ [2^{28} + 1, 2^{30}] & 2^{15} \\\\ [2^{30} + 1, 2^{32} - 1] & 2^{16} \\\\ \end{array} \end{aligned}

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 $\sqrt{x}$.

The Integer Square Root of the Integer Square Root

\begin{aligned} \left \lfloor \sqrt[4]{\alpha} \right \rfloor = \left \lfloor \sqrt{\left \lfloor \sqrt{\alpha} \right \rfloor} \right \rfloor \end{aligned}

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 $\lfloor \sqrt{\alpha} \rfloor = p$. Then, $p^2 \leq \alpha < (p + 1)^2$ since $p$ is an integer. Also, let $\lfloor \sqrt{p} \rfloor = q$. Then, $q^2 \leq p < (q + 1)^2$. So, $p + 1 \leq (q + 1)^2$. Therefore, $q^4 \leq p^2 \leq \alpha < (p + 1)^2 \leq (q + 1)^4$. That is, $q^4 \leq \alpha < (q + 1)^4$, which means $q \leq \sqrt[4]{\alpha} < q + 1$. Accordingly, $\left \lfloor \sqrt[4]{\alpha} \right \rfloor = q = \lfloor \sqrt{p} \rfloor = \left \lfloor \sqrt{\left \lfloor \sqrt{\alpha} \right \rfloor} \right \rfloor$.

Reference

[1] Henry S. Warren. 2012. Hacker’s Delight (2nd. ed.). Addison-Wesley Professional.