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:

1. Make an initial guess. Guess any positive number x0.
2. Improve the guess. Apply the formula x1 = (x0 + S / x0) / 2. The number x1 is a better approximation to sqrt(S).
3. 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.

1. An initial guess is x0 = 10.
2. Apply the formula: x1 = (10 + 20/10)/2 = 6. The number 6 is a better approximation to sqrt(20).
3. 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.

### How 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.

1. Count the integer digits in S. For example, S = 3249 has four digits.
2. 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.
3. 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.

Tags
Share Distinguished Researcher in Computational Statistics

Rick Wicklin, PhD, is a distinguished researcher in computational statistics at SAS and is a principal developer of PROC IML and SAS/IML Studio. His areas of expertise include computational statistics, simulation, statistical graphics, and modern methods in statistical data analysis. Rick is author of the books Statistical Programming with SAS/IML Software and Simulating Data with SAS.

1. Anders Sköllermo on

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.

2. 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 ?

```proc iml; / * begin IML session * / start MySqrt(x); / * begin module * / y = 1; / * initialize y * / do until(w<1e-3); / * begin DO loop * / z = y; / * set z=y * / y = 0.5#(z+x/z); / * estimate square root * / w = abs(y-z); / * compute change in estimate * / end; / * end DO loop * / return(y); / * return approximation * / finish; / * end module * /   t = MySqrt({3,4,7,9}); / * call function MySqrt * / s = sqrt({3,4,7,9}); / * compare with true values * / diff = t - s; / * compute differences * / print t s diff; / * print matrices * /```
• Thanks for the SAS/IML function. I didn't use SAS/IML because this operation does not involve vectors or matrices.

3. 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!

4. ryan derose on

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