  # Square root algorithm (example on while-loops)

 In this short article we’ll explore a square root algorithm as an excuse to use while-loops in our numerical software. We’re not going to use the built-in function 'sqrt'. Calculators typically implement routines to compute the exponential function and the natural logarithm, and then compute the root of a positive real number x using this identity: We’ll follow a different approach: an iterative method using while-loops.

The most common iterative method of square root calculation is known as the ‘Heron's method’ or ‘Babylonian method’.

First step, estimate a number. Think of this number as your first approach to a root (the closer to the actual square root of x, the fewer iterations will be needed to achieve the desired precision).

Second step, divide your number by this (known inexact) square root.

Third step, take the average of the result of the second step and ‘the root’.

Fourth step, use the result of the third step to repeat steps 2 and 3 until you have a number that is accurate enough for you.

The technique is a variation of the Newton-Raphson iterative solution method and it involves a simple algorithm, which results in a number closer to the actual square root each time it is repeated.

A code that accomplishes the algorithm is:

function ta = sr(x)
% Choose (arbitrarily) a first approach
fa = x/2;
% Divide your number by that first approach
sa = x/fa;
% Get the mean between the two previous numbers
ta = mean([fa sa]);

% Repeat until you obtain a good enough root
while abs(x - ta^2) > 0.000000001
fa = ta;
sa = x/fa;

ta = mean([fa sa])
end

We can use or try our function like this, from the command window for example:

format long
y1 = sr(10)
y2 = sqrt(10)

ta =    3.17857142857143
ta =    3.16231942215088
ta =    3.16227766044414
ta =    3.16227766016838

Our algorithm needed four iterations to obtain the final result (in y1, below). Note that only two iterations were needed to have an accuracy within three decimals. y2 is the result from the function ‘sqrt’, to compare with.

y1 =    3.16227766016838
y2 =    3.16227766016838

Another example,

y3 = sr(255)
y4 = sqrt(255)

ta =   34.34411196911197
ta =   20.88448273516376
ta =   16.54725251948965
ta =   15.97883290053478
ta =   15.96872262323155
ta =   15.96871942267163

Our algorithm needed six iterations to obtain the final result (in y3, below). Note that only four iterations were needed to have an accuracy within one decimal. y4 is the result from the function ‘sqrt’, to compare with.

y3 =   15.96871942267163
y4 =   15.96871942267131 