Square Roots in Python part 1: Estimates

Python provides us with a perfectly satisfactory method of calculating square roots, the math.sqrt() method, but in this pair of articles I will explore some of the ways square roots can actually be calculated from scratch. In this first article I will look into a few ways of calculating rough estimates which form a starting point or "seed" for the methods I'll examine in Part 2.

Estimating Square Roots

Many methods of calculating square roots require a starting value which is iteratively refined to home in on the final result. You can use a fixed arbitrary starting value but the more inaccurate it is the more iterations are required to get to an acceptably accurate final result.

It is therefore beneficial to calculate a very rough estimate to use as a starting value, and in this first article I will code a few ways of doing this.

The Project

This project consists of the following files:

• squarerootestimates.py
• squarerootplot.py
• squarerootestimatesdemo.py

The files can be downloaded as a zip, or you can clone/download the Github repository if you prefer.

Source Code Links

ZIP File
GitHub

Note that the ZIP file and Github repository also contain the files for Part 2 of this brace of articles.

The first file, squarerootestimates.py, implements a few different methods of calculating estimates.

In squarerootplot.py we plot the actual square roots of a range of numbers calculated by Numpy's sqrt method, together with any other values provided. These can be estimates, calculated values or both, and the plot gives us a quick and easy way to assess the accuracy of any values we may come up with.

Finally in squarerootestimatesdemo.py we call the functions in the previous two files to see which estimate function is best.

Using NumPy and Matplotlib

In this project I will be using NumPy instead of ordinary Python lists as it makes the calculations much more straightforward. If you aren't familiar with NumPy it basically provides us with a fast and efficient multi-dimensional numeric array and many associated methods. I'll only be using very basic NumPy functionality here which should be easy to understand. These are the links to NumPy's own site, and also its PyPI page which describes the PIP installation process.

https://numpy.org/

https://pypi.org/project/numpy/

I am also using Matplotlib for the plot, again using very basic functionality: we just throw a few values at Matplotlib and it will plot them for us in a nice little window. Another pair of links for you:

https://matplotlib.org/

https://pypi.org/project/matplotlib/

The Code

This is squarerootestimates.py which implements five methods of estimating square roots.

squarerootestimates.py

import numpy as np


def _nparray_to_sci_not(radicands):

exponents = np.log10(radicands).astype(int)

significands = radicands / 10**exponents

return {"significands": significands, "exponents": exponents}


def decimal_1(radicands):

sci_not = _nparray_to_sci_not(radicands)

estimates = 5.5 * 10**sci_not["exponents"]

return estimates


def decimal_2(radicands):

sci_not = _nparray_to_sci_not(radicands)

estimates = 3.16 * 10**sci_not["exponents"]

return estimates


def decimal_3(radicands):

sci_not = _nparray_to_sci_not(radicands)

filter_low = sci_not["significands"] < 10
    exponents_low = sci_not["exponents"][filter_low]

filter_high = sci_not["significands"] >= 10
    exponents_high = sci_not["exponents"][filter_high]

estimates_low = 2 * 10**exponents_low
    estimates_high = 6 * 10**exponents_high

return np.concatenate((estimates_low, estimates_high))


def linear_1(radicands):

sci_not = _nparray_to_sci_not(radicands)

estimates = (sci_not["significands"] - 1 + 1.2) * 10**sci_not["exponents"]

return estimates


def linear_2(radicands):

sci_not = _nparray_to_sci_not(radicands)

filter_low = sci_not["significands"] < 10
    significands_low = sci_not["significands"][filter_low]
    exponents_low = sci_not["exponents"][filter_low]
    estimates_low = (significands_low * 0.28 + 0.89) * 10**exponents_low

filter_high = sci_not["significands"] >= 10
    significands_high = sci_not["significands"][filter_high]
    exponents_high = sci_not["exponents"][filter_high]
    estimates_high = (significands_high * 0.89 + 2.8) * 10**exponents_high

return np.concatenate((estimates_low, estimates_high))

_nparray_to_sci_not

All the methods I have used here need the radicands (numbers for which we are calculating square roots) in scientific notation, for example 100 in scientific notation is 1 * 102. The first component (1 in this example) is called the significand and the second component (2 in this example) is called the exponent. In the _nparray_to_sci_not function we create two NumPy arrays, one for exponents and one for significands. These are then returned in a dictionary.

This function gives our first glimpse of the power and simplicity of NumPy: we can call functions and perform operations on an array as a single entity without a for/in loop.

decimal_1

This function takes a NumPy array of radicands (values to estimate the square roots of) and returns another NumPy array containing those estimates.

Firstly we call _nparray_to_sci_not to get the radicands scientific notation components. The exponents effectively break down the radicands into ranges: 1-10, 11-100, 101-1000 etc. We then multiply the radicands by 5.5 to get the mid-point as a rough estimate of the square root.

Again NumPy allows us to just use a single operation without a loop. I won't mention this again!

decimal_2

This is basically the same as the previous method but using 3.16 instead of 5.5

decimal_3

This is a refinement of the first two methods but using different multipliers depending on whether the radicand is < 10 or >= 10.

The radicands array is split into two parts, < 10 and >= 10, these are multiplied by 2 and 6 respectively, and then concatenated into a single result.

linear_1

This is the first of a more sophisticated estimate which provides a continuously variable set of estimates instead of a set of values fixed for each range of radicands.

Although it's more sophisticated it isn't particularly accurate so we'll move on to a better refinement.

linear_2

This method filters the significands and exponents of the radicands into two ranges, < 10 and >= 10. Estimates for each range are calculated by multiplying by and adding different values before concatenating the results into a single array.

Now let's move on to the code to plot square roots.

squarerootplot.py

import numpy as np
import matplotlib
import matplotlib.pyplot as plt


def plot(valuerange, sqrts=None):

matplotlib.pyplot.figure(figsize=(10,8), facecolor="#FFFFFF")

plt.title('Square Roots in Python ')

values = np.array(range(valuerange[0], valuerange[1]+1))
    npsqrts = np.sqrt(values)

plt.plot(values,
             npsqrts,
             label='NumPy sqrt',
             color='#0000FF',
             alpha=0.25,
             linewidth=4)

if sqrts != None:
        for sqrt in sqrts:
            plt.plot(values,
                     sqrt["sqrts"],
                     label=sqrt["label"],
                     color=sqrt["color"],
                     linewidth=1.0)

plt.legend()
    plt.show()

This function takes a range of radicands and an optional list of square roots which will be plotted. These actually consist of dictionaries containg a Numpy array of values and also a label and color.

Firstly we create a new Matplotlib figure of a certain size and color, and then set its title.

We then create a Numpy array of the required range, and calculate its square roots using Numpy's own sqrt method. These are plotted as definitive accurate square roots against which we can compare estimates or square roots calculated using other methods.

Next the list of dictionaries (if any) is iterated and the values plotted using the specified color and with the specified label.

Finally we call legend to show a key of labels, and show to actually display the plot.

Now we just need some code to call our estimate functions and plot the results.

squarerootestimatesdemo.py

import math

import numpy as np

import squarerootestimates as sre
import squarerootplot as srp


def main():

print("----------------------")
    print("| codedrome.com      |")
    print("| Square Roots       |")
    print("| Part  1: Estimates |")
    print("----------------------\n")

min = 1
    max = 99

estimates = []

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

estimates_decimal_1 = sre.decimal_1(radicands)
    estimates.append({"sqrts": estimates_decimal_1,
                      "label": "decimal_1",
                      "color": "#FF0000"})

estimates_decimal_2 = sre.decimal_2(radicands)
    estimates.append({"sqrts": estimates_decimal_2,
                      "label": "decimal_2",
                      "color": "#FF8000"})

estimates_decimal_3 = sre.decimal_3(radicands)
    estimates.append({"sqrts": estimates_decimal_3,
                      "label": "decimal_3",
                      "color": "#FF00FF"})

estimates_linear_1 = sre.linear_1(radicands)
    estimates.append({"sqrts": estimates_linear_1,
                      "label": "linear_1",
                      "color": "#00FFFF"})

estimates_linear_2 = sre.linear_2(radicands)
    estimates.append({"sqrts": estimates_linear_2,
                      "label": "linear_2",
                      "color": "#000080"})

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


if __name__ == "__main__":

main()

The min and max variable specify the range of radicands, estimates is a list of dictionaries to pass to squarerootplot.plot, and radicands is a Numpy array initialised to the specified range.

The five methods in squarerootestimates are then called. Each time a new dictionary is added to the estimates list, containing the estimates themselves as well as a label and color.

The squarerootplot.plot function is then called to plot the various estimates as well as the definitive square roots.

Now we can go ahead and run the program.

Running the Code

python3 squarerootestimatesdemo.py

If you haven't used Matplotlib before, it displays plots in a new window which you can resize and move around, and which also provides a simple toolbar at the bottom.

As I mentioned earlier the decimal estimates provide a fixed value for each range of radicands whereas the linear estimates vary continuously.

Obviously none of the estimates are accurate enough to be used as they are, but they provide a good starting point for the methods I'll implement in Part 2 of this series.

For updates on the latest posts please follow CodeDrome on Twitter