Integer Square Root

Newton’s Method

Newton

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αf(x) = x^2 - \alpha to know α\sqrt{\alpha} since it is a root of f(x)f(x). As a short review, Newton’s method computes iteratively to find the approximated root as follows. fn=fnxnxn+1    xnxn+1=fnfn    xn+1=xnfnfn    xn+1=xnxn2α2xn=xn2+α2xn=12(xn+αxn)\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 xnx_n is accurate to nn bits, then xn+1x_{n + 1} is accurate to 2n2n 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α23210 \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 xn+1=12(xn+αxn)x_{n + 1} = \lfloor \frac{1}{2} ( x_n + \lfloor \frac{\alpha}{x_n} \rfloor ) \rfloor. Also, xnx_n and α\alpha are positive integers. Then, the following statements hold.

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

Since xnx_n is an integer, the following holds. xn+1=12(xn+αxn)=12xn+αxn=12(xn+αxn)=xn2+α2xn\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 xn=(1+ϵ)αx_n = (1 + \epsilon) \sqrt{\alpha}. Noting that xn>αx_n > \sqrt{\alpha} since xn>αx_n > \lfloor \sqrt{\alpha} \rfloor and xnx_n is an integer, xn2+α2xn=xn+1xn2+α2xn    (1+ϵ)2α+α2(1+ϵ)α=xn+1<xn2+xn22xn    ϵ2+2ϵ+22(1+ϵ)α=xn+1<xn    2ϵ+22(1+ϵ)αxn+1<xn    αxn+1<xn\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 α1<xnα\sqrt{\alpha} - 1 < x_n \leq \sqrt{\alpha} since xn=αx_n = \lfloor \sqrt{\alpha} \rfloor. So, noting that xn2α<(xn+1)2x_n^2 \leq \alpha < (x_n + 1)^2, xn2+xn22xnxn2+α2xn=xn+1xn2+(xn+1)22xn    xnxn+1xn+1+12xn    αxn+1xn+1+12xn    αxn+1xn+1=α+1\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 xnx_n is larger than the true value α\lfloor \sqrt{\alpha} \rfloor, then the next guess xn+1x_{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 x0x_0 is equal to the least power of 2 that is greater than or equal to α\sqrt{\alpha}. For example, for α=4\alpha = 4, x0=2x_0 = 2, and for α=5\alpha = 5, x0=4x_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. Range of xInitial Guess[2,22]2[22+1,24]22[24+1,26]23[26+1,28]24[228+1,230]215[230+1,2321]216\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][0, 15], the expression (x > 0) + (x > 3) + (x > 8) computes x\sqrt{x}.

The Integer Square Root of the Integer Square Root

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


© 2024. All rights reserved.