Computing the Integer Square Root

Today I’m going to talk about a fast algorithm to compute the integer square root of a non-negative integer \(n\), \(\isqrt(n) = \lfloor \sqrt{n} \rfloor\), or in words, the greatest integer whose square is less than or equal to \(n\).^[1] Most sources that describe the algorithm take it for granted that it is correct and fast. This is far from obvious! So I will prove both correctness and speed below.

One simple fact is that \(\isqrt(n) \le n/2\), so a straightforward algorithm is just to test every non-negative integer up to \(n/2\). This takes \(O(n)\) arithmetic operations, which is bad since it’s exponential in the size of the input. That is, letting \(\Bits(n)\) be the number of bits required to store \(n\) and letting \(\lg n\) be the base-\(2\) logarithm of \(n\), \(\Bits(n) = O(\lg n)\), and thus this algorithm takes \(O(2^{\Bits(n)})\) arithmetic operations.

We can do better by doing binary search; start with the range \([0, n/2]\) and adjust it based on comparing the square of an integer in the middle of the range to \(n\). This takes \(O(\lg n) = O(\Bits(n))\) arithmetic operations.

// isqrt returns the greatest number x such that x^2 <= n. The type of
// n must behave like BigInteger (e.g.,
// https://github.com/akalin/jsbn ), and n must be non-negative.
//
//
// Example (open up the JS console on this page and type):
//
//   isqrt(new BigInteger("64")).toString()
function isqrt(n) {
  var s = n.signum();
  if (s < 0) {
    throw new Error('negative radicand');
  }
  if (s == 0) {
    return n;
  }

  // x = 2^ceil(Bits(n)/2)
  var x = n.constructor.ONE.shiftLeft(Math.ceil(n.bitLength()/2));
  while (true) {
    // y = floor((x + floor(n/x))/2)
    var y = x.add(n.divide(x)).shiftRight(1);
    if (y.compareTo(x) >= 0) {
      return x;
    }
    x = y;
  }
}

The core of the algorithm is the iteration rule: \[ x_{i+1} = \left\lfloor \frac{x_i + \lfloor \frac{n}{x_i} \rfloor}{2} \right\rfloor \] where the floor functions are there only because we’re using integer division. Define an integer-valued function \(f(x)\) for the right side. Using basic properties of the floor function, you can show that you can remove the inner floor: \[ f(x) = \left\lfloor \frac{1}{2} (x + n/x) \right\rfloor \] which makes it a bit easier to analyze. Also, the properties of \(f(x)\) are closely related to its equivalent real-valued function: \[ g(x) = \frac{1}{2} (x + n/x)\text{.} \]

For starters, again using basic properties of the floor function, you can show that \(f(x) \le g(x)\), and for any integer \(m\), \(m \le f(x)\) if and only if \(m \le g(x)\).

Finally, let’s give a name to our desired output: let \(s = \isqrt(n) = \lfloor \sqrt{n} \rfloor\).^[4]

(Note that any number greater than \(s\), say \(n\) or \(\lceil n/2 \rceil\), can be chosen for our initial guess without affecting correctness. However, the expression above is necessary to guarantee performance. Another possibility is \(2^{\lceil \lceil \lg n \rceil / 2 \rceil}\), which has the advantage that if \(n\) is an even power of \(2\), then \(x_0\) is immediately set to \(\sqrt{n}\). However, this is usually not worth the cost of checking that \(n\) is a power of \(2\), as is required to compute \(\lceil \lg n \rceil\).)

We are done proving correctness, but you might wonder if the check \(x_{i+1} \ge x_i\) in step \(3.2\) is necessary. That is, can it be weakened to the check \(x_{i+1} = x_i\)? The answer is “no”; to see that, let \(k = n - s^2\). Since \(n \lt (s+1)^2\), \(k \lt 2s + 1\). On the other hand, consider the inequality \(f(x_i) \gt x_i\). Since that would cause the algorithm to terminate and return \(x_i\), that implies that \(x_i = s\). Therefore, that inequality is equivalent to \(f(s) \gt s\), which is equivalent to \(f(s) \ge s + 1\), which is equivalent to \(g(s) = (s + n/s) / 2 \ge s + 1\). Rearranging yields \(n \ge s^2 + 2s\). Substituting in \(n = s^2 + k\), we get \(s^2 + k \ge s^2 + 2s\), which is equivalent to \(k \ge 2s\). But since \(k \lt 2s + 1\), that forces \(k\) to equal \(2s\). That is the maximum value \(k\) can be, so therefore \(n\) must be one less than a perfect square. Indeed, for such numbers, weakening the check would cause the algorithm to oscillate between \(s\) and \(s + 1\). For example, \(n = 99\) would yield the sequence \(\{ 16, 11, 10, 9, 10, 9, \dotsc \}\).

We will show that \(\NewtonSqrt\) takes \(O(\lg \lg n)\) arithmetic operations. Since each loop iteration does only a fixed number of arithmetic operations (with the division of \(n\) by \(x\) being the most expensive), it suffices to show that our algorithm performs \(O(\lg \lg n)\) loop iterations.

It is well known that Newton’s method converges quadratically sufficiently close to a simple root. We can’t actually use this result directly, since it’s not clear that the convergence properties of Newton’s method are preserved when using integer operations, but we can do something similar.

Define \(\Err(x) = x^2/n - 1\) and let \(ϵ_i = \Err(x_i)\). Intuitively, \(\Err(x)\) is a conveniently-scaled measure of the error of \(x\): it is less than \(1\) for most of the values we care about and it bounded below for integers greater than our target \(s\). Also, we will show that the \(ϵ_i\) shrink quadratically. These facts will then let us show our bound for the iteration count.

It’s also useful to show that if the initial guess \(x_0\) is bad, then the run-time degrades to \(Θ(\lg n)\). We’ll do this by defining the function \(\NewtonSqrt\) except that it takes a function \(\mathrm{I{\small NITIAL}\text{-}G{\small UESS}}\) that is called with \(n\) and assigned to \(x_0\) in step 1. Then, we can treat \(ϵ_0\) as a function of \(n\) and analyze how long \(ϵ_i\) stays above \(1\) to show that \(\NewtonSqrt\) uses an initial guess such that \(ϵ_0(n) = Θ(1)\), then Theorem 4 reduces to Theorem 3 in that case. However, if \(x_0\) is chosen to be \(Θ(n)\), e.g. the initial guess is just \(n\) or \(n/k\) for some \(k\), then \(ϵ_0(n)\) will also be \(Θ(n)\), and so the run time will degrade to \(Θ(\lg n)\). So having a good initial guess is important for the performance of \(\NewtonSqrt\)!