Python and successive approximation

I was doing some work in the yard and I wanted to know the smallest circle that would fit around a 4x6 inch post. Of course, just as a 2x4 is not 2 inches by 4 inches, a 4x6 post (what they call its "nominal" dimensions) is actually 3.5 inches by 5.5 inches, so this smallest circle has as its diameter the diagonal of that 3.5x5.5-inch rectangle. Pythagoras tells us that will be the square root of the sum of 3.5 squared and 5.5 squared. At the time, I only had my cell phone, which has a calculator with only basic math operations, so how do you get the square root of 42.5? Successive approximation.

You make your initial guess, knowing that it is greater than 6 but less than 7, and try 6.2 maybe, then 6.3, so on. The number 6.5 is pretty close with a square of 42.25, so you go to the next decimal place. Finally, you get as far as 6.5192, and say you're close enough. You don't need a square root function on your calculator for this.

Some time ago it occurred to me that the long-hand way of "calculating" a square root is nothing more than this same method on paper, so I thought it would be interesting to teach Python how to do this. Before I go any further, I would add that I am fully aware that the Python math module has a sqrt function for this, but that wouldn't be nearly as much fun as creating my own script.

My approach was to have the script ask for the value for which the square root is desired, then ask for the number of decimal places, and then go through the process of progressively approximating the square root, which is going to involve a lot of repetition, testing, and retesting. What I wanted then was some sort of incrementing loop, with a test to see if I was below or above the desired value. Here is the final result:

#!/usr/bin/env python
# sqrt.py - find square root

def sqtest(num, orig, ind):
    sq = 0
    while sq <= orig:
        num = num + ind
        sq = num * num
    num = num - ind
    return num

n = float(0) # this is our working number. We start at zero and increment.

innum = input("For what number would you like to get the square root? ")
i = input("To how many decimal places? ")

inc = float(10)
i += 1
if (innum < 0):
    i = 0
while i:
    inc = inc * 0.1
    n  = sqtest (n, innum, inc)
    i -= 1
if (innum >= 0):
    sq = n * n
    print "Approximate square root of " + str(innum)+ " is " + str(n) + "\n(Actual square of this would be "+str(sq)+")"
else:
    print "This would be an imaginary number"

For now, skip over the indented section, called a function, which begins with def sqtest. You're going to tell Python to begin the quest at zero, because you might at some point want the square root of a number between 0 and 1. After this, you ask for the number for which you need the square root, and then how many decimal places of precision you want.

Although it looks like your initial increment is going to be 10, this will be refactored so you can smoothly work with decimals later. Add 1 to the number of decimal places because you will have a cycle for units and above before you get to the decimals. Next, you check to see if the number input was less than zero, so you don't waste time on imaginary numbers.

Now comes the loop, where you send your working number n, your original number innum, and the current increment, initially 10 x 0.1, off to the sqtest function.

In sqtest, to which you send the values for your current working number, the number for which you are needing a square root, and the current incrementing value. After clearing any prior value, you check the square of the working number plus the new increment, and if it's lower than innum, increment until the square is greater than innum, at which point you back up one incremental notch, and send back the new working number. The while loop lower down now takes your last increment value and takes one-tenth of that, and you're off to sqtest again. You repeat until you've finished the requested number of decimal places.

Here is a sample output:

For what number would you like to get the square root? 28837465
To how many decimal places? 6
Approximate square root of 28837465 is 5370.052606
(Actual square of this would be 28837464.9912)

Notice that I not only display the square root, but also show what the actual square of that would be, for a better sense of how close I am. One thing to keep in mind is the uncertainty of the last digit. Because of the way this script works, this value is truncated at 6 decimal places, and not rounded, so you don't know if the seventh decimal place would be greater or less than 5. It's better to have more places than you think you need, so that any rounding you might decide to do is accurate.

You might think that starting from 0 and incrementing by 1 would be time-consuming for large numbers, but I think you'll be surprised how fast the interpreter runs here—the square root of this 8-digit number with 6 decimal places was pretty instantaneous. For real-world numbers, this is quite adequate, and even with 15-digit numbers plus 6 decimal places, it took perhaps 3 seconds. However, at some point, there might be issues with the number of significant digits.