Python and successive approximation

Learn how to do successive approximation in Python with this example script.
398 readers like this.
Best of Opensource.com: Programming

Zagrev on Flickr. CC BY-SA 2.0

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.

Greg Pittman
Greg is a retired neurologist in Louisville, Kentucky, with a long-standing interest in computers and programming, beginning with Fortran IV in the 1960s. When Linux and open source software came along, it kindled a commitment to learning more, and eventually contributing. He is a member of the Scribus Team.

7 Comments

Here is some homework. With just a few modifications, you can turn this into a cube root finder.
Here is the new output:

For what number would you like to get the cube root? 384857
To how many decimal places? 6
Approximate cube root of 384857 is 72.738855
(Actual cube of this would be 384856.992155)

It's called the bisection method.

Bisection sounds like cutting something in half, but in reality it's successive approximation.

In reply to by Zahir Jacobs (not verified)

One of my thoughts as a neurologist is that another real-world kind of successive approximation would be a pitcher practicing pitching a ball. This is why, if you have the skills (code?) you get progressively better with practice.

I do something similar with a ksh (Korn shell) script called dec2dyadic. It converts a decimal real number to an integer plus a dyadic fraction (ie, the denominator is a power of two) using a similar sort of iteration. It's good for being able to take a real number (20.777) and turn it into something on a standard ruler or yardstick (20 25/32). Or, you can have it take some number (length), and similarly split it up into an integral number of dyadic parts (like cutting a board into 5 equal lengths). And you can specify a margin or sensitivity. This was before I started learning python. ;-) You can find a copy here: ftp://ftp.cs.duke.edu/pub/des/scripts/dec2dyadic

I'm just learning Python right now and enjoyed this. I've made the adjustments for it to work with Python3:

#!/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 = eval(input("For what number would you like to get the square root? "))
i = eval(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")

Creative Commons LicenseThis work is licensed under a Creative Commons Attribution-Share Alike 4.0 International License.