# The Root of All Squares

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: n^{2}=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.

Here's a bit of output from a sample run:

```
Square root of 3368 = 58 58.034187 Check=3367.966797 error=-0.033203 (-0.000986%)
Square root of 7739 = 88 87.971428 Check=7738.972168 error=-0.027832 (-0.000360%)
Square root of 12 = 3 3.428571 Check=11.755102 error=-0.244898 (-2.040815%)
Square root of 6226 = 79 78.904457 Check=6225.913574 error=-0.086426 (-0.001388%)
Square root of 8586 = 93 92.659462 Check=8585.776367 error=-0.223633 (-0.002605%)
Square root of 8094 = 90 89.966484 Check=8093.968262 error=-0.031738 (-0.000392%)
Square root of 7539 = 87 86.826591 Check=7538.856934 error=-0.143066 (-0.001898%)
```

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

/* Old fashioned square root calculation -- Williams */ /* Compile: gcc -D TEST -o sqrt1 sqrt1.c */ /* This relies on the series: n^2==1+3+7+...+(2n-1) The ENIAC computed square roots this way. Other than test code this requires no multiplication and only division by 2. */ /* Define TEST if you want a TEST program Define ROOT_FLOAT if you want a floating point answer passed back through a pointer to float */ #include <stdio.h> #ifdef TEST #include <stdlib.h> /* random numbers */ /* If testing, must have float! */ #define ROOT_FLOAT #endif #ifdef ROOT_FLOAT int eniac_isqrt(int n0, float *v) #else int eniac_isqrt(int n0) #endif { int n=n0; /* original number */ int odd; /* the odd number we are subtracting */ int adjust; /* which way to adjust (depends on which is closest) */ int err1, err2; /* error terms */ #ifdef ROOT_FLOAT float interpolate; /* fractional linear interpolation */ #endif /* First subtract odd numbers until we hit zero or negative */ for (odd=1;;odd+=2) { n-=odd; if (n<=0) { if (n==0) adjust=1; /* odd is exactly 2n-1 */ break; } } /* n is the amount we went negative, if any */ if (n!=0) /* not an exact 2n-1 so work out which is closer */ { err1=abs(n+odd); err2=abs(n); adjust=err1>err2?1:-1; } #ifdef ROOT_FLOAT /* find the fractional part through linear interpolation */ if (v && n!=0) { if (err1>err2) { interpolate=-err2; } else { interpolate=err1; } interpolate/=odd; if (v) *v=(odd+adjust)/2.0+interpolate; } else if (v) *v=(odd+adjust)/2.0; /* or if exact... */ #endif /* return integer result */ return (odd+adjust)>>1; /* note >>1 is the same as /2 */ } #ifdef TEST /* Simple routine to drive the calculation and print stats */ void debugsqrt(int x) { float v; int iv=eniac_isqrt(x,&v); printf("Square root of %d = %d %f",x,iv,v); #ifdef TEST v*=v; printf(" Check=%f error=%f (%f%%)\n",v,v-x,((v-x)/x)*100); #else putchar('\n'); #endif } /* Simple main */ // generate random numbers from 0 to MAXTESTINT-1 #define MAXTESTINT (10000) // Random number seed change for different sequence of randoms */ #define SEED 1 void main(int argc, char *argv[]) { int i; srand(SEED); debugsqrt(9); /* known test case ! */ for (i=0;i<100;i++) { debugsqrt(rand()%MAXTESTINT); } } #endif