Square Roots in Python part 2: Calculations

In the first part of this article I wrote code to calculate rough estimates of square roots which enable formulas which converge on accurate square root values to do so with fewer iterations. Part 1 also included code to graph estimates to compare them to accurate values.

In this second part I'll implement some of those formulas, again graphing the results alongside definitive values.

Calculating Square Roots

There are many methods of calculating square roots and for this post I've implemented a representative sample of three of the better known ones: Babylonian, Bakhshali and Exponential Identity. I'll describe each one before we get to coding them.

The Babylonian Method

This is believed to be the oldest method, and is also known as Heron's Method after Hero of Alexandria. We start with an estimate (as calculated in part 1) and then divide the radicand (number we want the square root of) by the estimate, and take the average of the two values. This forms the next estimate and is more accurate than the previous one. Specifically, if the first estimate is too high the new one will be too low and vice versa. We then repeat the process until we get an acceptably accurate result.

The following is an example of four iterations of the Babylonian Method to find the square root of 100, which of course is 10. The estimates are bold so you can watch them getting ever closer to 10. Note that the values alternate between being too high and too low, as mentioned in the previous paragraph.

The Babylonian Method for Radicand S = 100
Estimate 1 (Linear Method)	36
S / Estimate 1	2.77777777777778
Estimate 2 = Mean of previous 2 values	19.3888888888889
S / Estimate 2	5.15759312320917
Estimate 3 = Mean of previous 2 values	12.273241006049
S / Estimate 3	8.14780708296315
Estimate 4 = Mean of previous 2 values	10.2105240445061
S / Estimate 4	9.79381661157797
Estimate 5 = Mean of previous 2 values	10.002170328042

The Bakhshali Method

This method is a refinement of the Babylonian method and although it's more complex and uses two intermediate values, a and b, it does converge quicker than the Babylonian so requires fewer iterations.

This is the Bakhshali Method in action, again with a radicand of 100.

The Bakhshali Method for Radicand S = 100
Estimate 1 (Linear Method)	36
a = (S - Estimate12) / (2 * Estimate1)	-16.6111111111111
b = Estimate1 + a	19.3888888888889
Estimate2 = b - (a2 / (2 * b))	12.273241006049
a = (S - Estimate22) / (2 * Estimate2)	-2.06271696154294
b = Estimate2 + a	10.2105240445061
Estimate3 = b - (a2 / (2 * b))	10.002170328042

Exponential Identity

This method is by far the simplest to implement and needs no initial estimate or iterations, although it does require the calculation of a logarithm which is abstrated away by a Python function.

This is it in mathematical notation.

Clearly we can implement this with just a single line of Python.

The Project

This project consists of the following files:

• squarerootestimates.py
• squarerootplot.py
• squarerootcalculations.py
• squarerootcalculationsdemo.py

The files can be downloaded as a zip, or you can clone/download the Github repository if you prefer. The ZIP and repository contain the files for both parts of this article.

Source Code Links

ZIP File
GitHub

The first two files, squarerootestimates.py and squarerootplot.py, are listed and described in Part 1.

The squarerootcalculations.py file contains functions to calculate square roots and squarerootcalculationsdemo.py shows the previous files in use.

NumPy and Matplotlib

This project uses Numpy arrays instead of Python lists as this streamlines the code. I have also used Matplotlib to plot the graphs.

The Code

This is squarerootcalculations.py.

squarerootcalculations.py

import math

import numpy as np

import squarerootestimates as sre


def babylonian_array(radicands, iterations):

estimates = sre.linear_2(radicands)

vfunc = np.vectorize(babylonian)
    square_roots = vfunc(radicands, estimates, iterations)

return square_roots


def bakhshali_array(radicands, iterations):

estimates = sre.linear_2(radicands)

vfunc = np.vectorize(bakhshali)
    square_roots = vfunc(radicands, estimates, iterations)

return square_roots


def exponential_identity_array(radicands):

vfunc = np.vectorize(exponential_identity)
    square_roots = vfunc(radicands)

return square_roots


def babylonian(s, estimate, iterations):

sqrt = estimate

for i in range(0, iterations):

sqrt = (sqrt + s / sqrt) / 2

return sqrt


def bakhshali(s, estimate, iterations):

sqrt = estimate

a = 0
    b = 0

for i in range(0, iterations):

a = (s - sqrt**2) / (2 * sqrt)

b = sqrt + a

sqrt = b - (a**2 / (2 * b))

return sqrt


def exponential_identity(n):

return math.e ** (0.5 * math.log(n))

babylonian_array, bakhshali_array and exponential_identity_array

These three functions take a NumPy array of radicands (values to calculate the square root of) and return another NumPy array of the square roots.

The actual square roots are calculated by separate functions which take individual radicands. We therefore need to use NumPy's vectorize method which creates a wrapper function for the specified function; this wrapper function can operate on and return NumPy arrays.

The first two methods of calculating square roots, Babylonian and Bakhshali, need estimates so these are calculated using the linear_2 method from the squarerootestimates module.

babylonian, bakhshali and exponential_identity

These calculate the square root of the radicand s, the first two using the specified estimates and number of iterations. They are straightforward implementations of the methods described earlier.

squarerootcalculationsdemo.py

import math

import numpy as np

import squarerootcalculations as src
import squarerootplot as srp


def main():

print("------------------------")
    print("| codedrome.com        |")
    print("| Square Roots         |")
    print("| Part 2: Calculations |")
    print("------------------------\n")

babylonian()

# bakhshali()

# exponential_identity()


def babylonian():

min = 1
    max = 99

plot_data = []

radicands = np.array(range(min, max + 1))

square_roots_babylonian = src.babylonian_array(radicands, 0)
    plot_data.append({"sqrts": square_roots_babylonian,
                      "label": "Babylonian 0 iterations",
                      "color": "#000000"})

square_roots_babylonian = src.babylonian_array(radicands, 1)
    plot_data.append({"sqrts": square_roots_babylonian,
                      "label": "Babylonian 1 iterations",
                      "color": "#FF0000"})

square_roots_babylonian = src.babylonian_array(radicands, 2)
    plot_data.append({"sqrts": square_roots_babylonian,
                      "label": "Babylonian 2 iterations",
                      "color": "#0000FF"})

square_roots_babylonian = src.babylonian_array(radicands, 3)
    plot_data.append({"sqrts": square_roots_babylonian,
                      "label": "Babylonian 3 iterations",
                      "color": "#00FF00"})

square_roots_babylonian = src.babylonian_array(radicands, 4)
    plot_data.append({"sqrts": square_roots_babylonian,
                      "label": "Babylonian 4 iterations",
                      "color": "#FF8000"})

srp.plot((min,max), plot_data)


def bakhshali():

min = 1
    max = 99

plot_data = []

radicands = np.array(range(min, max + 1))

square_roots_bakhshali = src.bakhshali_array(radicands, 0)
    plot_data.append({"sqrts": square_roots_bakhshali,
                      "label": "Bakhshali 0 iterations",
                      "color": "#000000"})

square_roots_bakhshali = src.bakhshali_array(radicands, 1)
    plot_data.append({"sqrts": square_roots_bakhshali,
                      "label": "Bakhshali 1 iterations",
                      "color": "#FF0000"})

square_roots_bakhshali = src.bakhshali_array(radicands, 2)
    plot_data.append({"sqrts": square_roots_bakhshali,
                      "label": "Bakhshali 2 iterations",
                      "color": "#FF8000"})

srp.plot((min,max), plot_data)


def exponential_identity():

min = 1
    max = 99

plot_data = []

radicands = np.array(range(min, max + 1))

square_roots_exponential_identity = src.exponential_identity_array(radicands)
    plot_data.append({"sqrts": square_roots_exponential_identity,
                      "label": "Exponential Identity",
                      "color": "#FF8000"})

srp.plot((min,max), plot_data)


if __name__ == "__main__":

main()

The main function simply calls three other functions which we can comment/uncomment to run one at a time.

The babylonian function creates a NumPy array of radicands between 1 and 99. It then calls babylonian_array five times with ever-increasing iterations, including 0 iterations which will just return the estimate provided. Each time the result is added to the plot_data list. Each result is in a dictionary which also contains a label and color.

I have used fixed numbers of iterations here, the number required to give a reasonably accurate square root. If I were coding a real life implementation I'd square the result after every iteration and check it against the radicand using a tolerance value, say 0.01, and terminate the iterations if it passed the test.

Finally we call squarerootplot.plot to show the values after various numbers of calculations.

The bakhshali and exponential_identity functions work in the same way, but note that I have only gone up to 2 iterations with Bakhshali, and the Exponential Identity method is a one-hit calculation that does not require iterations.

Now we can run the code:

Running the Code

python3 squarerootcalculationsdemo.py

This is the graph for the Babylonian method.

The topmost, and therefore least accurate, line is black and is the result of passing 0 as the iterations argument, thereby just getting back the estimates provided. As you can see the red, blue, green and orange lines get ever closer to the thick light blue line which shows the definitive square root values as calculated by NumPy. The final orange line matches the NumPy values to a high degree of accuracy.

The green 3-iterations line is very close to the light blue and orange lines and at this scale seems to touch them so let's zoom in on a part of the plot so we can distinguish them more easily.

This shows just the part of the plot for radicands 6 to 16, and here we can see that the green line does actually diverge from the light blue and orange lines.

In main comment out babylonian(), uncomment bakhshali() and run the program again.

The Bakhshali method has converged to a very accurate result after just two iterations, but remember the calculations are more complex so this does not necessarily mean the function is more efficient.

Now do the comment/uncomment thing once more to run the exponential_identity function.

Of course there are no estimates or iterations here, just a single set of values which fit with the NumPy light blue values very accurately. Bear in mind though that this method does need to calculate the logarithm of the radicand so is more than a simple piece of arithmetic.

For updates on the latest posts please follow CodeDrome on Twitter