HOWTO-square-root (2001/08/15)

I thought it was about time to publish a nifty, little algorithm for computing the square-root of a number. Here it is:

• Split the number in pieces of two digits each (ex.: 65536 becomes 6'55'36)
• From the leftmost digit(s), start subtracting the odd numbers beginning with 1 and count how many times you could do it w/o underflow (6: -1 -3, two times)
• The number of times you could do it is the first digit of the square-root. Multiply it by 10, add the successor and combine the remainder of the previous subtraction with the next two digits (2 * 10 = 20; 20 + 21 = 41; remainder was 6 - 4 = 2, next two digits 55, operate on 255)
• Now subtract the odd numbers beginning with the abovementioned sum and again count how often you can do it (255: -41 -43 -45 -47 -49; 5 times and the remainder is 255 - (5*45) = 30)
• Again multiply the so-far-result by 10, add to its successor and append the next two digits of the original number for the next subtraction loop (25*10 = 250; 250 + 251 = 501; remainder was 30, next two digits 36, operate on 3036)
• And once more subtract the sequence of odd numbers beginning with the sum of (partial result) * 10 + (partial result) * 10 + 1. (3036: -501 -503 -505 -507 -509 -511; 6 times and there is no remainder)
• The square-root of 65536 is 256

Do you somehow doubt this algorithm actually works and do you think the lines above are a hoax? Ok, then it's time for another proof that you can verify with your calculator ;-)

• Finding the square-root for 1234
• Split each two digits: 12'34
• Start subtracting odd numbers from the leftmost digit-packet 12 without underflowing: -1 -3 -5 works 3 times, remainder is 12 - 9 = 3
• 3 * 10 = 30; 30 + 31 = 61; remainder (3) with next two digits appended: 334
• 334: -61 -63 -65 -67 -69; worked 5 times, remainder 334 - 325 = 9
• 35 * 10 = 350; 350 + 351 = 701
• Oh! We now have no more two-digit-packets left, but we have a remainder!?
  Simple: we set the decimal point in the result now.
  Result so far was 35 so we have 35, and now take double zeroes to extend the remainder and continue.
  Remainder was 9, append 00, so operate on 900 now.
• 900: -701; works once, remainder 900 - 701 = 199 (solution so far 35,1)
• 351 * 10 = 3510; 3510 + 3511 = 7021
• 19900: -7021 -7023; works twice, remainder 19900 - 14044 = 5856 (solution so far 35,12)
• 3512 * 10 = 35120; 35120 + 35121 = 70241
• 585600: -70241 -70243 -70245 -70237 -70249 -70251 -70253 -70255; worked 8 times, remainder 585600 - 562000 = 23600
• I'll stop here, because the numbers are too long for my poor little brain to work on ;-)

The square-root of 1234 is approx. 35,128.... Is that accurate enough for a quick guess?

You can run this algorithm without a calculator, at least to some degree, even if you are a math illiterate (like me ;-)

Now for the best part of it all: This algo works just as well if implemented for base 16. Instead of two digits you take two nibbles (also known as a byte) and instead of multiplying by 10 you multiply by 16 (or left shift four times). This does sound like you could implement it pretty damn fast in any assembly language, doesn't it? No multiplication, only simple subtraction loops...

I wanted to contribute this algorithm to the GnuMP project, but I have yet to find the time and devotion to dig into the GMP number format to implement it correctly. And I would also have to think about when to stop the calculation of irrational square-roots. The algorithm lets you continue until infinity, if you like, but you probably don't have that much time :-)

Continued...

I thought a little bit longer about the algorithm and found a way to compute the cubic root of a number in a similiar manner. Take a look at the tables below. It seems you can find a cubic root by subtracting the sequence of step * 3rd odd numbers (1, 7, 19, 37 etc.) from each triplet of digits? But I don't (yet) know what to do with the remainders, ie. what sum do I have to build to begin the subtraction with in the next step:

• 11*11*11 = 1331 -> 1'331
• 1 -1; works once, remainder 0, next digits 331; result so far 1, but how do I get from here to 331 to start the next step with?
• 331 -331 = 0; works once, result 11, remainder 0
• The cubic root of 1331 is 11 .. I know it, but how do I do it? :-)

12*12*12 = 1728; remainder 728 which is 331 + (331+66)
13*13*13 = 2197; remainder 1197 which is 331 + (331+66) + (331+66+72)
14*14*14 = 2744; remainder 1744 which is 331 + (331+66) + (331+66+72) + (331+66+72+78)
15*15*15 = 3375; remainder 2375 which is 331 + (331+66) + (331+66+72) + (331+66+72+78) + (331+66+72+78+84)
etc.

Now what is so special about the magic number 331? And where do the additional summands for the following values come from? What's the origin of this number magic and why is the additional distance in every new term 66 + 6 * (step-1)? I don't see the background for this.. only the facts.

It looks like I have to go back to the drawing board (or the keyboard for that matter) and play around with the numbers to see if and how a cubic root can really be found with an algorithm similiar to the one for the square root.

Table of square, cubic and "to the power of 4" numbers

There is an order in the 3rd column of distances too, do you see it? Look at the distances between the distances of the n-th odd numbers... 7-1=6, 25-7=18 (6+12), 55-25=30 (6+12+12), 97-55=42 (6+12+12+12), 151-97=54 (6+12+12+12+12), 217-151=66 (6+5*12) etc. I wonder how this could be implemented in an algorithm... and I also wonder if one could generalize this approach to computing roots of any order?

More information

After asking on the usenet group sci.math.num-analysis I was pointed in the right direction and was told, that the algorithm was known and used in mechanical calculators ("Friede desk calculator") in the early 60s at Boeing already. Later Hewlett Packard sold calculators such as the HP9100, which used a similiar or identical algorithm utilizing BCD (binary coded decimals).

I found this reference to the inner workings of the HP Model 9100A and the mathematical background for the algorithm. The background is that:

   n
   Σ  2*i - 1 = n2
  i=1

so it says that the sum of all 2*i-1 for i starting with 1 up to n is equal to the square of n. I would be interested to learn how to write the corresponding formula for the cubic numbers, any hints?

It also seems someone or some people thought it was not possible to transport this algorithm to the binary world. But as I pointed out above, it is not only possible, it is outright simple to do it. So IMHO there is no good reason not to use this algorithm in math libraries, because it should be a lot faster than using Newton's approximation (which involves at least one division for every approximation step), or even more complicated approaches based on logarithms - I guess (IANAM).

Have fun

Juergen