## DSP Trick: Square Root Computation

Date: Mon, 19 Apr 1999 17:50:23 +0200
From: rainer storn <rainer.storn@hl.siemens.de>
Newsgroups: comp.dsp
Subject: Square root computation
THIS WORK IS PLACED IN THE PUBLIC DOMAIN

**Name:** Square root computation

**Category:** Algorithmic Trick

**Application:** e.g. Magnitude computation

**Advantages:** Simple and fast

**Introduction:**

Once in a while the square root of data has to be computed on a DSP. A
method that comes to mind is polynomial approximation of the square root
function. A much simpler method, however is the usage of Newton's method or, if
division is a problem, successive approximation.

**The Trick:**

1) Newton's method:

Newton's method states that if you want to find the zero of a function f(x)
you may start at some guessing point x1 and compute the next guessing point for
the zero x2 according to

x2 = x1 - f(x1)/f´(x1) (1)

In the next step make x1 := x2 and repeat equation (1) etc. until
the difference between x1 and x2 is sufficiently small. Newton's method does
not alway converge, but if it does, it converges fast. In the special case of
the square root computation Newton's method is guaranteed to converge.

Let's say you want to compute x = sqrt(z) then make

f(x) = x*x - z. (2)

If we compute f(x) = 0 we immediately see, that x must be the
square root of z.

Now plug (2) into (1) to get

x2 = x1 - 0.5*x1 - 0.5*z/x1
= 0.5*(x1 - z/x1);

which is the recursive equation to be used.

The starting value for x1 is not critical but the closer it is to the real
zero the better. I usually use x1 = z.

2) Successive approximation:

If division is a problem, then successive approximation might be the right
choice to compute sqrt(x). Successive approximation basically does a binary
search to find the square root. A little C-style pseudocode might best explain
how it is done:

//----y = sqrt(y)---------
a = 2^ceil(b/2); //b is the wordlength you are using
y = a; //a is the first estimation value for sqrt(x)
for (i=0; i<ceil(b/2); i++) //each iteration achieves one bit
{ //of accuracy
a = a*0.5;
diff = x - y*y;
if (diff > 0)
{
y = y + a;
}
else if (diff < 0)
{
y = y - a;
}
}

Again, there might be better choices for the starting point but the one
given is safe.

Rainer Storn