Flashback: I found the Friden algo part of this article printed and abandoned (wasted) in a common room in uni. Took home one copy (There were few for some reason) and loved it. Today I found embedded.com article on HN & looked for this article. And thought HNers would love.

This is excellent advice.

See http://groups.google.fr/group/comp.lang.c/browse_thread/thre...

The method is very simple. We know how to compute the square of a value. So start by finding the biggest integer value with a single bit set whose square is smaller than the input value. Then set the bit just below to one, compute the square of the value and compare it with the input value. If it is bigger, the bit should be 0 otherwise it should be left to 1. And so on.

The thing is to avoid the multiplication to compute the square. So if you have a value n and you want to set a bit to one. We compute in fact the square of (N+B) where B is an integer with the bit to set to one, set to one. We then have (N+B)² = N² + 2NB + B² The last term is easy to compute from the previous B with the previous bit set. We simply divide by 4. 2NB is computed by a simple shift and N² is the result of the previous computation.

For fixed point, we benefit from the fact that while we progress in computing the square root, we have an increasing number of bits set to 0 in the most significant bits of the values. I then simply shift up to get usable bits in the less significant bits.

The only limitation is that we need the sign bit of the fixed point. So the computation won't work with big unsigned fixed points. This is solved by unwinding the first iteration.

It took me some time to understand all this. All I had was a short code to compute integer square root and I had no idea how it worked.

In the vicinity of the correct root, the method converges exponentially, meaning that it gives twice as many digits of accuracy at each step

should instead be converges quadratically

In fact it's somewhat interesting, that quadratic error decay equals exponential precision growth:

  0.1^2 = 0.01
  0.01^2 = 0.0001
  0.0001^2 = 0.00000001
  ...

"As an iterative method, the order of convergence is equal to the number of terms used. With 2 terms, it is identical to the Babylonian method; With 3 terms, each iteration takes almost as many operations as the Bakhshali approximation, but converges more slowly. Therefore, this is not a particularly efficient way of calculation."

Heh, woops. My second inkling was the way slide rules do it: exponentials and logarithms ( sqrt(n)=exp(.5*ln(n)) ) Turns out that's more on target:

"Pocket calculators typically implement good routines to compute the exponential function and the natural logarithm, and then compute the square root of S"

Yeah, they are accurate, but don't need to be fast. The blink of an eye is a long time.