Compute square root using (bit) additions and shifts as primitives

Question

Question: Given an $n$-bit natural number $N$, how to compute $\lceil \sqrt{N} \rceil$ using only $O(n)$ (bit) additions and shifts?

The tip is to use binary search. However, I could not achieve the required complexity (I got $O(n^2)$).

---

What does it mean by using only $O(n)$ (bit) additions and shifts:

This is an exercise in an algorithm book.
In my opinion, it means that adding two, say $n$-bit, natural numbers costs $O(1)$ and shifting a, say $n$-bit, natural number also costs $O(1)$. Then we are only allowed to use such $O(1)$ operations $O(n)$ times.
It does not mention the cost of comparison. I guess we can ignore it or assume that comparing two, say $n$-bit, natural numbers costs $O(1)$ as well.

---

My $O(n^2)$ algorithm:

1. Determine the range of the number of bits $t$ of $\lceil \sqrt{N} \rceil$: $$2^{\frac{n-1}{2}} \le \sqrt{N} \le 2^{\frac{n}{2}} \Rightarrow 2^{\lfloor \frac{n-1}{2} \rfloor} \le \lceil \sqrt{N} \rceil \le 2^{\lceil \frac{n}{2} \rceil}$$ Therefore, $$t_1 \triangleq \lfloor \frac{n-1}{2} \rfloor + 1 \le t \le \lceil \frac{n}{2} \rceil + 1 \triangleq t_2.$$
2. Binary search: Find $\lceil \sqrt{N} \rceil$ between $2^{t_1}$ and $2^{t_2}$ using binary search. For each number $x$, to compute $x^2$ using additions and shifts as primitives and compare it with $N$.

The complexity is thus $O(n \times n) = O(n^2)$ for $O(n)$ times of binary search and computing $x^2$, each of which in turn takes $O(n)$ additions and shifts.

Answer 1

An iterative algorithm seems like it should work.

Let $M=\lfloor N/4 \rfloor$. Suppose we know that $x$ is the integer approximation to $\sqrt{M}$, i.e., $x=\lceil \sqrt{M} \rceil$, and suppose we know the value of $x^2$ (obtained previously).

Now we want to find $y=\lceil \sqrt{N} \rceil$. What are the possible values of $y$? I'm pretty sure the only possible values are $y=2x$ or $y=2x+1$. And, it's easy to try both of them and see which is correct. In particular, for $y=2x$, we have $y^2=4x^2$, which can be obtained from $x^2$ by two left-shifts ($O(1)$ time); for $y=2x+1$, we have $y^2=4x^2+4x+1$, which can be obtained from $x^2$ and $x$ with four left-shifts and two additions ($O(1)$ time). Now just compare those two values to $N$ to see which one is correct.

In this way, we get an iterative algorithm where we do $n/2$ iterations, and where each iteration takes $O(1)$ time. The total running time is $O(n)$, as required.

I realize this didn't use binary search. Oh well.

Answer 2

Are we talking integers here? Where N is n bits long?

A = 2(n/2), B = A  and C = A2
Step: B = B/2
     If C > N,  
         C = C - 2AB + B2    // too high - make smaller
         A = A - B
     Else 
         C = C + 2AB + B2   // keep this bit
         A = A + B                 
Repeat until B = 0                  // =1 on last loop

Loop is performed n/2 times, which should give you O(n) performance

Edit: How does it work, & why?
This is a version of Successive Approximation, which is also used in CORDIC algorithms.
Starting with the largest possible single bit (with a square less than N) you set one bit at a time, and calculate the new square.
If the new square is still less than N, keep the bit as set.
If the new square is too big, clear the bit, undo the effect of adding it, and move on to the next bit.

Example: N = 441 (1 1011 1001 binary), n = 9

Start:  A = 24 = 16 (1 0000)  B = 16 C = 256 (100 0000)

1   B = 8 (1000) C = 256 + 2(16)(8) + (8)(8) = 576 (10 0100 0000) {high}
    A = 16 + 8 = 24
2   B = 4  (100) C = 576 - 2(16)(4) + (4)(4) = 400 (1 1001 0000) {low}
    A = 24 - 4 = 20
3   B = 2   (10) C = 400 + 2(20)(2) + (2)(2) = 484  (1 1110 0100) {high}
    A = 20 + 2 = 22
4   B = 1    (1) C = 484 - 2(20)(1) + (1)(1) = 441  (1 1011 1001) {keep this}
    A = 22 - 1 = 21
5   B = 1/2 or 0 in integer math; end

Answer 3

The main method is to fill in the bits of $\sqrt{N}$ from left to right while keeping our estimate below it, or rather the square of our estimate below $N$. Each bit $b$ is a power of 2, so squaring or multiplying another number by $b$ is always a bit shift.

If the current estimate is $a$, $b=2^i$, and we know $a^2$ already, we get $(a+b)^2 = a^2 + 2ab + b^2$, and we can rewrite the second and third terms as $ a << (i+1)$ and $1 << (i << 1)$. We then add it all up and test (I assume you can do $<$) and set bit $i$ if the square is still below $N$.

We start the loop at $i = n/2 = n >> 1$ and count down to zero, keeping $a$ and $a^2$ as we go. It's a kind of binary search, but one where the bounds map to single-bit differences.

Answer 4

I like Alan Campbell's answer: with careful tracking of previous guesses, the new subtraction is easy each time, and the binary shift-and-add square root is about as fast as a binary shift-and-add divide.

But it may be possible to go faster, by instead of making your next guess a single binary digit, instead using an "A-b" x "A-b" algorithm, and making your next guess the average of your previous guess, and the original number divided by the previous guess. That sounds like it would take longer, not shorter. However, the division doesn't have to be exact. So if the division only runs to the square root of the number of digits remaining to find, then you might actually save time. Moreover, if for your division you use the french method, of shorthand division, then you might actually break some speed in your calculation for really big divides.

Now, if we add in calculations in parallel that yield preliminary correctable results before the answer is found... then we might be onto something.