When I was in the sixth grade, I learned an iterative procedure for computing square roots by hand. Yes, I said by hand. Scientific calculators with a square root key were not yet widely available, so I and previous generations of children suffered through learning to calculate square roots by hand.
I still remember being amazed when I first saw the iterative square root algorithm. It was the first time that I thought math is magical. The algorithm required that I make an initial guess for the square root. I then applied a "magic formula" a few times. The magic formula improved my guess and estimated the square root that I sought.
The Babylonian square-root algorithm
The iterative method is called the Babylonian method for finding square roots, or sometimes Hero's method. It was known to the ancient Babylonians (1500 BC) and Greeks (100 AD) long before Newton invented his general procedure.
Here's how it works. Suppose you are given any positive number S. To find the square root of S, do the following:
- Make an initial guess. Guess any positive number x0.
- Improve the guess. Apply the formula x1 = (x0 + S / x0) / 2. The number x1 is a better approximation to sqrt(S).
- Iterate until convergence. Apply the formula xn+1 = (xn + S / xn) / 2 until the process converges. Convergence is achieved when the digits of xn+1 and xn agree to as many decimal places as you desire.
Let's use this algorithm to compute the square root of S = 20 to at least two decimal places.
- An initial guess is x0 = 10.
- Apply the formula: x1 = (10 + 20/10)/2 = 6. The number 6 is a better approximation to sqrt(20).
- Apply the formula again to obtain x2 = (6 + 20/6)/2 = 4.66667. The next iterations are x3 = 4.47619 and x4 = 4.47214.
Because x3 and x4 agree to two decimal places, the algorithm ends after four iterations. An estimate for sqrt(20) is 4.47214.
The Babylonians and Greeks could estimate square roots by hand. Algorithm and #SAS program. Click To TweetHow do you choose an initial guess?
You can choose any positive value as an initial guess. However, when I was a kid I used to race my friends to see who could find the square root the fastest. I discovered that if you choose a good guess, then you have to compute only a few iterations. I invented a strategy for finding a good guess, which I call "The Rule of Twos and Sevens."
The Rule of Twos and Sevens chooses an initial guess from among the candidates {2, 7, 20, 70, 200, 700, ...}. The Rule use the fact that a square root has about half as many integer digits as the number itself. The Rule follows and is illustrated by using S = 3249.
- Count the integer digits in S. For example, S = 3249 has four digits.
- Consider the candidates that have half as many digits as S. (If S has an odd number of digits, round up the number of digits.) For example, 20 and 70 are two-digit candidates. Use those candidates when the number S has three or four integer digits.
- Use mental arithmetic to decide which candidate is better. For example, 20 is the square root of 400 whereas 70 is the square root of (about) 5000. Because S is closer to 5000, choose 70 as the starting guess. (Challenge: Convince yourself that you can use 20 for all three-digit integers, 200 for all five-digit integers, and so forth.)
In other words, a good guess starts with a 2 or a 7 and has about half as many digits as are in the whole-number part of S. My experience is that The Rule of Twos and Sevens usually converges to a solution (within two decimal places) in four or fewer iterations.
For small numbers, you can choose the initial guess more precisely. If S is less than 225 (=15x15), you can bracket the square root by using the perfect squares {1, 4, 9, ..., 196, 225} and then use the corresponding integer in the range [1, 15] as an initial guess.
Implement the Babylonian square-root algorithm in SAS
In tribute to all students who ever struggled to perform this algorithm by hand, I implemented the Babylonian square-root algorithm in SAS. You can implement the algorithm directly as a DATA step program, but I chose to use PROC FCMP to define new DATA step functions.
proc fcmp outlib=work.funcs.MathFuncs; /* Rule of 2s and 7s: Count the number of digits in S. Choose initial guess to have half as many digits (rounded up) and start with a 2 or 7. */ function BabylonianGuess(S); /* provide initial guess */ str = put(floor(S), 16.); /* convert [S] to string */ L = length(strip(str)); /* count how many digits */ d = ceil(L/2); /* about half as many digits (round up) */ guess2 = 2*10**(d-1); guess7 = 7*10**(d-1); if abs(guess2**2 - S) < abs(guess7**2 - S) then return( guess2 ); else return( guess7 ); endsub; /* the Babylonian method (aka, Hero's method) for finding square roots */ function BabylonianSqrt(S); epsilon = 100*constant("maceps"); /* convergence criterion */ if S < 0 then x = .; /* handle negative numbers */ else if S=0 then x = 0; /* handle zero */ else do; x = BabylonianGuess(S); /* initial guess */ xPrev = 0; do while (abs(x - xPrev) > epsilon); xPrev = x; x = (xPrev + S/xPrev)/2; /* iterate to improve guess */ end; end; return( x ); endsub; quit; /* Functions defined. Compare Babylonian algorithm to modern SQRT function */ options cmplib=work.funcs; /* define location of functions */ data Compare; input S @@; BabySqrt = BabylonianSqrt(S); /* Babylonian algorithm */ Sqrt = sqrt(S); /* modern computation */ Diff = abs(BabySqrt - Sqrt); /* compare values */ datalines; -3 0 1 2 4 10 16 30 100 3249 125348 ; proc print label; label BabySqrt = "Babylonian Sqrt"; run; |
Notice that the Babylonian method produces essentially the same numerical value as the built-in SQRT function in SAS. Of course, the iterative BabylonianSqrt function is not nearly as efficient as the built-in SQRT function.
For more information about the Babylonian algorithm and other algorithms for computing square roots, see Brown (1999) Coll. Math J.
Do you remember learning a square root algorithm in school? Do you still remember it? Are there other "ancient" algorithms that you have learned that are now unnecessary because of modern technology? Leave a comment.
14 Comments
Hi! Third root solution: X(n+1)= (X(n) + S/X(n)**2) *1/2 which gives X**3= S.
The main question is: Did the Babylonians know this and did they use it.
Rick,
Why not post a IML function ? Actually in IML documentation there is already such a function. I thought it was writen by yourself wasn't it ?
Thanks for the SAS/IML function. I didn't use SAS/IML because this operation does not involve vectors or matrices.
Great post Rick, makes me think how lucky we are to have all this computing power at our fingertips; after all, Computer used to be a job title!
Pingback: All I really need to know about Newton's method I learned in primary school - The DO Loop
Pingback: Compute the square root matrix - The DO Loop
Pingback: The math you learned in school: Yes, it’s useful! - The DO Loop
I think you would get better guesses in your system by using 6 instead of 7.
An even-digit number lies in the range [10^(n-1), 10^n), e.g. when a number has 4 digits, it is at least 1,000 and less than 10,000.
The geometric mean of these extremes is 10^((2n - 1)/2),
The square root of which is 10^((2n-1)/4) = 10^(n/2) * 10^(-1/4)
When n is odd, you get a power of 10 times 1.778... , which is close to 2
but when n is even, you get a power of 10 times 5.623... , which is much closer to 6 than to 7
8943 this 4 digits number become 0097. 8289. 0618. 5954.. 5619. 7909 how is calculated
(8943+483 ) using babilonian method
Please answer me
My method for finding the initial guess is simply to find the closest perfect square and then estimating how far your number is from that perfect square. I mean, even from 1 to 10 this covers alot..... 1 4, 9, 16, 25, 36, 49, 64, 81, 100,
It's kind of easy for 4 digits because you just count by 10s, so 100s when squared 100, 400, 900 1600, etc
If you wanted to do a number like 550 or something then you are going somewhere between 400 and 900. (answer between 20 and 30 but closer to 20)
I mean I'm really overexplaining here because with your number , 3249, it doesn't take a lot of effort to see that 3600 is the closest perfect square to that number.
Since 3249 is smaller than 3600, you know your square root is going be smaller than 60. You can also somewhat tell that it is not a whole smaller than 60, and in fact if go down to the previous, of 2500, and try to guage where 3249 sits between 2500 and 3600, you would just apply a similar ratio between 50 and 60.
Half way between 2500 and 3600 would be about 3050 so you know 3249 is just over half way between 50 and 60.
I mean you can probably just look at 3249 and not think too hard to see that 56 is probably a fairly reasonable initial guess, but go ahead and try it and tell me how many iterations you need.
Thanks for writing. Your method is known as "linear interpolation."
Do you know a litterature reference for the formula which generates the continues fraction for (1 + x)^0.5, please? There is a lot to be found in the web, regarding that fraction, but I have not been able to find out who came up with it... The formula that I have in mind is
(1 + x)^0.5 = 1 + x / (1 + (1 + x)^0.5) (1)
On the basis of that formula, I have used
(1 + x)^0.5 # (3x + 4) / (x + 4) (2)
after approximating
(1 + x)^0.5 # 1 + x/2 (3)
in the numerator in Eq.1. In other words, I contented myself with a single iteration of Eq.1.
Eq. 2 gives (1 + x)^0.5 with a 1% error or less, for 0<x<1, by the way, whereas the error in Eq.3 is 6% for x=1.
Thank you
No, I don't know the name of this equation. It is a tautology, because
sqrt(1+x) = 1 + x/(1 + sqrt(1+x))
I think you can derive it by using the continued fraction algorithm. From this equation, you can recursively replace the sqrt(1+x) term on the right side to go to any depth.
You might want to study the continued fraction algorithm by starting with the continued fraction for sqrt(2). See https://www.cut-the-knot.org/proofs/SqContinuedFraction.shtml.
See if you can replace '2' with '1+x' and perform the same steps.
I would think that rather than guessing between decimal multiples of 2 and 7, you would extract the even powers of 10 and then look for the square root of numbers between 1 and 100, (or extract all powers of ten, and if its odd (say 2n+1), you calculate the square root of a number between 1 and 10 and then multiply it by tabulated value of sqrt(10)* 10^n. Anyways, this was just a thought. One can probably improve the babylonian method in a variety of ways.