Calculating a Square Root

In Terry Pratchett’s book Hogfather, a series of hilarious events leads a professor named Ponder Stibbons to wait up for the Disc World’s version of Santa Claus.  “I’m waiting for the Hogfather, thought Ponder Stibbons.  I’m in the dark waiting for the Hogfather.  Me.  A believer in Natural Philosophy.  I can find the square root of 27.4 in my head. (He’d have to admit that the answer would be ‘five and a bit,’ but at least he could come up with it.)  I shouldn’t be doing this.”

This begs the question: how do you find the square root of any number in your head?  Like Ponder Stibbons, we can easily find the closest square number and get a good estimate: 27.4 is just over 25, which is 5², so the square root of 27.4 is “5 and a bit.”  But is there an algorithm for finding the square root of any number mentally, or with paper and a pencil?  These days it’s so easy to mindlessly plug in the problem into a nearby calculator (on a cell phone, tablet, or computer) that we don’t stop to think about how it can be done.

There are actually several different ways of calculating a square root of a number, but we will focus on one particular method.  First, start by guessing an answer.  (The closer the guess is to the correct answer, the less steps we will need to take in the long run.)  Let’s guess that the square root of 27.4 is 5.  Second, we divide our guess into the number being square-rooted.  In our example, 27.4 divided by 5 is 5.48.  Therefore, we know that the square root of 27.4 is somewhere between 5 and 5.48.  Third, we will take the average of those two numbers (5.24) and repeat the process with 5.24 as our new guess.  Here is the result of a few rounds of calculations:

Guess	Divide	Average
5	27.4 / 5 = 5.48	5.24
5.24	27.4 / 5.24 = 5.229007634	5.234503817
5.234503817	27.4 / 5.234503817 = 5.234498046	5.234500931
5.234500931	27.4 / 5.234500931 = 5.234500931	5.234500931

In this case, after the second round of calculations we have an answer that is accurate to five decimal places, and after the third round of calculations we have an answer that is accurate to at least nine decimal places.

Some of us may be able to do the first round of calculations for finding a square root in our heads, but even Professor Ponder Stibbons would have a difficult time with subsequent rounds.  In our example above, dividing 27.4 by 5.24 is no easy task (although it can be done).  Now the question becomes: how can we program a computer to find a square root using only the basic mathematical functions?

Algebraically, we can express this process as an iterated function.  Let p equal the number to be square rooted, x equal the guess, and y equal the square root.  Then y is the average of ^p/_x and x, or y = (^p/_x + x) / 2, or y = ^{p + x^2}/_2x.  In our case above, p = 27.4 and x₀ = 5, which means x₁ = ^{27.4 + 5^2}/_2∙5 = 5.24, x₂ = ^{27.4 + 5.24^2}/_2∙5.24 = 5.234503817, and so on.  The actual square root would be the fixed point of this iterated function, which can be found by letting y = x.  Therefore, y = ^{p + y^2}/_2y, and multiplying both sides by 2y, 2y² = p + y², and subtracting y² from both sides, y² = p or y = √p.  Since we have expressed this process as an iterated function, the first round’s solution can be represented as x₁ = (p + x₀²)/(2x₀), the second round’s solution as x₂ = (p + x₁²)/(2x₁), the third round’s solution as x₃ = (p + x₂²)/(2x₂), and so on.  Using substitution and a little bit of algebra, we can deduce that x₂ = (p² + 6px₀² + x₀⁴)/(4x₀(p + x₀^2)) and x₃ = (p⁴ + 28p³x₀² + 70p²x₀⁴ + 28px₀⁶ + x₀⁸)/(4x₀(p + x₀²)(p² + 6px₀² + x₀⁴)), where each equation is a closer and closer approximation to y = √p.  In fact, using x₀ = 5, the graph of x₃ and √p are nearly indistinguishable from about 0 < p < 250, and x₃ > √p when p > 250.

Graph of x₃ when x₀ = 5 and y = √p

(0 < p < 900, 0 < y < 30)

The equation for x₂ or x₃ may give a good approximation of a square root for numbers less than 250, but we ought to have something that works well for all numbers.  Fortunately, iterated functions are loops, and pretty much all computer programs can loop the same set of commands until a certain criteria is met.  In this case, we would want the computer to loop through the formula ^{p + x^2}/_2x for subsequent terms of x until the same decimal approximation is reached for two subsequent terms.  The following is a program that does exactly that in the computer language Python:

01	#Square Root Function
02	#Python 2.7.3
03
04	#function to determine if a string is a number
05	def is_number (n):
06	#if the number can successfully change into a float then return true
07	try:
08	float(n)
09	return True
10	#otherwise return false.
11	except:
12	return False
13
14	#function to determine if a number is positive
15	def is_positive (n):
16	#first check to make sure n is really a number
17	if (is_number(n)):
18	#if the number is less than zero then return false
19	if (float(n) < 0):
20	return False
21	#otherwise return true
22	else:
23	return True
24	#if n is not a number return false
25	else:
26	return False
27
28	#function to find the square root without using a built-in square root function
29	def square_root (p):
30	#first check to make sure p is a positive number
31	if is_positive(p):
32	#assign variables
33	p = float(p)
34	r = 9 #number of decimal places to round to
35	x_old = 5 #starting guess for x0
36	#loop values in y = (p + x^2)/(2x) until x_old = x_new (rounded)
37	while True:
38	#find x_new using y = (p + x^2)/(2x)
39	x_new = (p + x_old * x_old) / (2 * x_old)
40	#round x_old and x_new
41	x_old_round = "{:.{prec}f}".format(x_old, prec=r)
42	x_new_round = "{:.{prec}f}".format(x_new, prec=r)
43	#if x_old and x_new are equal, then stop and return value
44	if (x_old_round == x_new_round):
45	return float(x_new_round)
46	break
47	#make x_old the x_new value and start the loop over again
48	x_old = x_new
49	#if p is not a positive number then return an error message
50	else:
51	return "ERROR"

Running the above program and typing in the command to use the square_root function for 27.4 results in the following answer:

>>	square_root(27.4)
5.234500931

This accurate and quick result only needed about 50 lines of code (and only 30 lines if you take out all the comments).

Therefore, there is an easy algorithm for finding the square root of a number.  If you are good at mental math and can remember the formula ^{p + x^2}/_2x, you could even impress your friends by finding the square root to the nearest hundredth place (well, at least Ponder Stibbons would be impressed).  But the best approximation for a square root requires several iterations of the formula ^{p + x^2}/_2x which can be done quickly with a short computer program.