Al Williams

Dr. Dobb's Bloggers

The Root of All Squares

May 06, 2011

I've written this algorithm in an obscure assembly language before, but I decided to write a version in C to demonstrate how simple it is to take a square root using this method.

There have been several times that I have had trouble stifling a chuckle while conducting a job interview. I remember one time I was interviewing a prospective embedded programmer and it became clear he was primarily a C programmer. I asked him how he would calculate a square root. He told me about the sqrt function!

I was actually looking for a discussion about Newton's approximation. I would have accepted some talk of CORDIC or some other algorithm too. But a call to a library wasn't what I had in mind (it was that kind of job). The project was on very anemic hardware and I eventually rewrote an old algorithm that only used addition, subtraction, and division by two. It is the same technique used by the old ENIAC computer, so if it works on such ancient hardware, it should be a good fit for a simple CPU.

The idea is to find the sum of the odd integers that is just greater than the number. So for the square root of 9, for example, 1+3+5+7 is the smallest sum of odds that is greater than 9. Taking the last number, 7, you can determine that the square root is between (7-1) /2 and (7+1) /2. In other words, between 3 and 4. Of course, the answer is 3. This works because of the formula: n2=1+3+5+....+(2n-1). Actually, if you recognize that 9=1+3+5 (that is, the result is zero, not negative) you can stop at the 5 and compute n=(5+1)/2.

This even works when the result is not an integer if you are willing to settle for an integer estimate or interpolate the answer. Take 30, for example. Subtracting, we find the first odd integer that throws us into negative territory is 11 (the left over is -6 for 11 and +5 for 9). So we know the answer is a little closer to 5 than it is to 6. The integer approximation is, therefore, 5 (the right answer is about 5.48). We can do better using linear interpolation. Taking the error value (5) and dividing by the number (11) we get 0.45 and adding that to the 5 gets us pretty darn close (5.45 versus 5.48).

There are even further optimizations possible, especially if you expect to handle larger numbers (see http://is.gd/KsFqLR for a good treatment of your options). Also, if you are using fixed-point math anyway, don't forget you can multiply the number by a power of 2 (by shifting left) and then divide the answer by the square root of the original power (which will also be a power of 2, and therefore just needs right shifts). This will increase precision without the need for interpolation. For instance, you might multiply the input number by 64 and then divide the result by 8.

As an example of using the powers of 2 to scale the input and output, imagine you want to find the square root of 129. If you multiply it first by 64 (6 left shifts), you get 8256. The integer square root algorithm will report the root is 91. Now divide 91 by 8 (the square root of 64) to get 11.375 (which squares to 129.391, about a 0.3% error). My calculator puts the real answer at about 11.358. The algorithm with interpolation yields 11.348.

I've written this algorithm in an obscure assembly language before, but I decided to write a version in C to demonstrate how simple it is to take a square root using this method. The actual subroutine only uses adding, subtracting, and right shifting for dividing by two. Porting it should be easy. If you want a floating-point result, define ROOT_FLOAT, but then you do get a floating point division (it only happens twice per square root; once to interpolate and once to divide by 2).

The test program (define TEST when compiling) assumes you are also defining ROOT_FLOAT and it takes more liberties, doing multiplications to check the subroutine, divisions to compute error, and calling the C library for random number generation. Of course, you don't have to port any of that to a new target.

The output shows the integer result followed by the floating-point result. Then it shows the result of squaring the floating point result, the raw error, and the error as a percentage. In general, smaller numbers that don't work out exactly are going to have larger errors (note the square root of 12 is only off by 0.245 but that's over 2%).

For the purists, keep in mind that the code computes the positive square root, but that making it negative results in an equally valid square root! And, of course, if you need an exact answer, you'd need a better algorithm or to use the result as a starting point and iterate through nearby values until you found one that met your error criteria.

The code appears below. Although I compiled it on Linux, there's nothing mysterious about it and it should run just about anywhere. After all, if the ENIAC could handle it....

Dr. Dobb's encourages readers to engage in spirited, healthy debate, including taking us to task.
However, Dr. Dobb's moderates all comments posted to our site, and reserves the right to modify or remove any content that it determines to be derogatory, offensive, inflammatory, vulgar, irrelevant/off-topic, racist or obvious marketing or spam. Dr. Dobb's further reserves the right to disable the profile of any commenter participating in said activities.

Video

This month's Dr. Dobb's Journal

This month,
Dr. Dobb's Journal is devoted to mobile programming. We introduce you to Apple's new Swift programming language, discuss the perils of being the third-most-popular mobile platform, revisit SQLite on Android
, and much more!