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