Since it sometimes is required to calculate the square root of some number. This job is straight forward by solving the 2nd order equation for the square root of 2 i.e. \[ \begin{equation} x^2 − 2 = 0 \Rightarrow x = \sqrt{2} \approx 1.41 \end{equation} \]

However, when we need to define how a computer is suppose to calculate the square root, we can use numerical methods. I finished my B.Sc.EE project this semester. The project theme was Embedded Signal Processing with the design and implementation of a GFSK-Demodulator for Satellite Communication, in a fixed point DSP. For the design it was necessary to calculate the Magnitude \[M\] of a complex number \(I + jQ\) which is defined as \(M = \sqrt{I^2 + Q^2}\). Therefore, I came across several interesting algorithms for approximation of the square root in fixed point arithmetic. I believe this is an obvious subject for my first blog post here, where I will demonstrate some of these in Python, starting with Newtons Method.

Newtons Method is an iterative approach, with an initial guess, \(s\), as seen here

\[ s = \frac{1}{2}(s + x/s)\]

where \(x\) is the number from which it is desired to calculate the square root.

As an example, if we want to calculate the square root of 25 we can make a loop like this.

In [15]:

x = 25.
s = 1.
kmax = 6
for k in range(kmax):
    s = 1./2. * (s + x/s)
    print "After %s iterations, s = %20.15f" % (k+1,s)

After 1 iterations, s =   13.000000000000000
After 2 iterations, s =    7.461538461538462
After 3 iterations, s =    5.406026962727994
After 4 iterations, s =    5.015247601944898
After 5 iterations, s =    5.000023178253949
After 6 iterations, s =    5.000000000053722

Newtons Methods approaches the correct answer rather quick, compared to the initial guess being \(1\), which is relatively far from \(5\). An automatic guess could be the fraction i.e

\[ s = \frac{x}{4}\]

In [16]:

x = 25.
s = x/4.
kmax = 6
for k in range(kmax):
    s = 1./2. * (s + x/s)
    print "After %s iterations, s = %20.15f" % (k+1,s)

After 1 iterations, s =    5.125000000000000
After 2 iterations, s =    5.001524390243903
After 3 iterations, s =    5.000000232305737
After 4 iterations, s =    5.000000000000005
After 5 iterations, s =    5.000000000000000
After 6 iterations, s =    5.000000000000000

Then the square root approximation algorithm has reached the correct \(5\) already after 4 iterations, which is quite fair.

Another situation is when we have a fixed point machine for implementing such an algorithm. In such a situation we have to work primarily with numbers which is a power of \(2\). Any mulitplication or division with a power of two is normally done without mulitplying, since it can be done by a shifting the binary registers from inside the ALU.

When calculating the complex magnitude \(M = \sqrt{I^2 + Q^2}\) as described above a similar procedure can be utilized. To approximate the square root operation, the Binary-Shift Magitude Estimation can be used to estimate \(M\).

\[ M \approx A\text{M}_{ax} + B\text{M}_{in}\]
where \(A=\frac{15}{16}, B = \frac{15}{32}\) which is similar sizes as used in Newtons Method

In [17]:

import math
I = 3.
J = 11.
A = 15/16.
B = 15/32.

if I > J:
    M = A*I + B*J
else:
    M = B*I + A*J
print "M = %10.5f" %M
M2 = math.sqrt(I**2 + J**2)
print "M = %10.5f" %M2

M =   11.71875
M =   11.40175

This algorithm’s appeal is that it requires no explicit multiples because the \(A\) and \(B\) values are binary fractions and the formula is implemented with simple binary right-shifts, additions, and subtractions. (For example, a 15/32 times \(I\) multiply can be performed by first shifting \(I\) right by four bits and subtracting that number from \(I\) to obtain 15/16 times z. That result is then shifted right one bit). This is simulated in the following script.

In [18]:

I = 3
J = 11
J4=J>>4; J5=J>>5; I4=I>>4; I5=I>>5;

if I > J:
    M = (I-I4) + ((J-J5)>>1)
else:
    M = (J-J4) + ((I-I5)>>1)
print "M = %10.5f" %M
## For comparison
M2 = math.sqrt(I**2 + J**2)
print "M2 = %10.5f" %M2

M =   12.00000
M2 =   11.40175

Many combinations of \(A\) and \(B\) are provided in scientific litterature, yielding various accuracies for estimating \(M\). All these algorithms are open to modification and experimentation. This will make them more accurate and computationally expensive, or less accurate and computationally less demanding. Your accuracy and data throughput needs determine the path you take to a root of less evil.

Further reading about this subject can be suggested here:

A Root Less Evil - DSP tips & tricks

Lastly, a small piece of square root code is submitted.

In [18]:

In []:

My Sqrt mysqrt.py download

'''
Module for approximating the square root
Input
        x: The number for which the square root is approximated.
output
        s:  The requested square root.
'''
def newton(x,debug=False):
    '''newton(x)'''
    
    from numpy import nan

    x = float(x)
    if x==0.:
        return 0.
    elif x<0.:
        print '*** Error, x must be nonnegative'
        return nan
    assert x>0. and type(x) is float,"Unrecognized input"
    
    s = x/5.
    kmax = 100
    tol = 1.e-14
    for k in range(kmax):
        if debug:
            print "Before iteration %s, s = %20.15f" % (k,s)
        s0 = s
        s = 1./2. * (s + x/s)
        delta_s = s - s0
        if abs(delta_s / x) < tol:
            break
    if debug:
        print "After %s iterations, s = %20.15f" % (k+1,s)
    return s

def test():
    from numpy import sqrt
    import mysqrt

    xvalues =[0., 2., 100., 10000, 1.e-4]
    for x in xvalues:
        print 'Testing with x = %20.15e' % x
        s = mysqrt.newton(x)
        s_numpy = sqrt(x)
        print 's = %20.15e,   numpy.sqrt = %20.15e' \
                % (s, s_numpy)
        assert abs(s-s_numpy) < 1.e-14, \
                'Disagree for x = %20.15e' % x

python4enlightening

the open notebook

A root less evil

Comments