Computing Square Roots

The Museum of HP Calculators

HP Forum Archive 19

Computing Square Roots
Message #1 Posted by Thomas Klemm on 27 Mar 2009, 5:43 p.m.

While this might be obvious to some of you I only recently realized that the method used to calculate the square root in probably most HP-calculators is in fact closely related to Horner's method.
Let's assume we want to calculate the square root of 41 thus solving the equation:

x² - 41 = 0.
First we apply a trick which can be found in Cochran's patent and multiply the whole equation by 5:
5x² + 0x - 205 = 0
Now let's consider the following simple coordinate transformations:

shift: x' = x - 1
zoom: x' = 10 * x
What happens to a quadratic equation if we apply these transformations?

Shift

P(x) = ax² + bx + c = Q₁(x)(x-1) + c' = (Q₂(x)(x-1) + b')(x-1) + c' = (a'(x-1) + b')(x-1) + c' = (a'x' + b')x' + c'
Thus by continued polynominal division we can get the coefficients of the transformed polynominal.
We can use Horner's method to calculate the remainders.

Zoom

ax² + bx + c = 0 | * 100 100ax² + 100bx + 100c = 0 ax'² + 10bx' + 100c = 0 a'x'² + b'x' + c' = 0
By comparing the coefficents we get:
a' = a b' = 10b c' = 100c

Example
Here I assume you're familiar with the Horner scheme:

5 0 -205
1 5 5 -200 5 10 5 -------------------------------------- 5 10 -200 step
1 5 15 -185 5 20 5 -------------------------------------- 5 20 -185 step
1 5 25 -160 5 30 5 -------------------------------------- 5 30 -160 step
1 5 35 -125 5 40 5 -------------------------------------- 5 40 -125 step
1 5 45 -80 5 50 5 -------------------------------------- 5 50 -80 step
1 5 55 -25 5 60 5 -------------------------------------- 5 60 -25 step
1 5 65 40 here we turn the wheel back
5 600 -2500 and zoom
From now on I've chosen a shorter format:
5 610 -1895 step 5 620 -1280 step 5 630 -655 step 5 640 -20 step 5 650 625 back 5 6400 -2000 zoom 5 6410 4405 back 5 64000 -200000 zoom 5 64010 -135995 step 5 64020 -71980 step 5 64030 -7955 step 5 64040 56080 back 5 640300 -795500 zoom 5 640310 -155195 step 5 640320 485120 back 5 6403100 -15519500 zoom 5 6403110 -9116395 step 5 6403120 -2713280 step 5 6403130 3689845 back 5 64031200 -271328000 zoom
Now Cochran's trick becomes evident: In the second column appears the result 6.40312.
We could have also used a separate counter but it's not necessary.

What are we actually doing?
We move our coordinate system step by step to the right until we step over the root. Then we go one step back and zoom. Thus our unit becomes one tenth of what it was before. Then we continue ...
In pseudocode:
for (1 .. 12): while c <= 0: step back zoom
The same algorithm could be used to solve the cubic root or x² = x + 1. Unfortunately Cochran's trick won't work in these cases.
Personally I consider this algorithm a beautiful pearl. It's quiet amazing that the same algorithm was used in all these calculators over the years with only minor changes.

Links

HP35 square root algorithm

J.E. MEGGIT, Pseudo Division and Pseudo Multiplication processes, IBM Journal, Res & Dev, April 1962, 210-226.

US Patent 3576983: DIGITAL CALCULATOR SYSTEM FOR COMPUTING SQUARE ROOTS


Re: Computing Square Roots
Message #2 Posted by Bob Patton on 27 Mar 2009, 9:00 p.m.,
in response to message #1 by Thomas Klemm

This post reminded me of something I wrote for my kids a while back. It mighty be of passing interest to this forum.
Square Roots
When I was going to school, one mark of an educated person was the ability to find a square root with pencil and paper. While personal calculators have largely eliminated that need, you might still have a problem finding an exact root of a very large number, say 30 digits.
I doubt that they teach the process any more. It's somewhat akin to long division. In fact, I wonder if it's not already essentially lost except for those who study ancient methods. I decided to try and document it. It was hard to remember and it's fairly difficult to describe in words alone. I came up with a fairly compact description that relies on a somewhat algebraic description. Here it is.
Square Root - pencil/paper digit by digit method
NTBR = Number to be rooted ASF = answer so far, ignoring any decimal point BAL = balance
Group NTBR in pairs from the decimal point.
Start with left-most pair (or single left-over digit.) Find largest n such that n^2 less-than-or-equal the pair. Subtract n^2 from the pair to get remainder. This n is the first digit of the answer.
Repeat for as many digits as wanted:
Create next BAL by appending next-2-digits to remainder. If next BAL is ever 0 you have an exact square root at that point.
Then find largest n such that (20ASF + n)n less-than-or-equal BAL.
Subtract (20ASF+n)n from BAL to get remainder. This n is the next digit of the answer.
You can set up the mechanics of this somewhat like the long division process with pencil and paper
Example: square root of 125
1 1. 1 8 0 3 ________________ \/1 25.00 00 00 00 1 1^2 = 1 -- 2_|0 25 1| 21 <-- 21 X 1 +---- 22_ | 4 00 1 | 2 21 <-- 221 X 1 +----- 222_ |1 79 00 8 |1 78 24 <-- 2228 X 8 +------- 2236_ | 76 00 0 | 0 <-- 22360 X 0 +-------- 22360_ | 76 00 00 3 | 67 08 09 <-- 223603 X 3 +--------- 8 91 91


Re: Computing Square Roots
Message #3 Posted by Thomas Klemm on 28 Mar 2009, 7:24 a.m.,
in response to message #2 by Bob Patton

No wonder this reminds you of the traditional digit-by-digit method since essentially they are the same. It's just a different way to look at it.
Your example using Horner scheme:
1 0 -125
10 1 10 -25 1 20 1 -------------------------------------- 1 20 -25
1 1 21 -4 1 22 1 -------------------------------------- 1 220 -400
1 1 221 -179 1 222 1 -------------------------------------- 1 2220 -17900
8 1 2228 -76 1 2236 1 -------------------------------------- 1 22360 -7600
0 1 22360 -7600 1 22360 1 -------------------------------------- 1 223600 -760000
3 1 223603 89191 1 223606 1


Re: Computing Square Roots
Message #4 Posted by Karl Schneider on 28 Mar 2009, 8:43 p.m.,
in response to message #1 by Thomas Klemm

Thomas and Bob --
Thank you for posting the traditional square-root technique and the links.
I was taught the "paired digit" technique in school, but we didn't spend much time on it, and all I remembered was the pairing of digits and the first step.
I didn't know that there was an IBM Journal. Perhaps it was the inspiration for Hewlett-Packard Journal.
-- KS

Re: Computing Square Roots
Message #5 Posted by Michael de Estrada on 28 Mar 2009, 9:33 p.m.,
in response to message #4 by Karl Schneider

I went to college in the 1960's pre-calculator days, and we were equipped with a slide rule and a book with tables of trig and log functions. The "trick" we were taught to find a root of a number was to use logarithms as follows:
nth root of X = anti-log [(log X)/n]
Based on the relationship log a^b = b * log a
Hence, the cube root of 27 = anti-log [(log 27)/3] = 3
I do remember in elementary school learning long division, however, I can't remember learning to find square roots. Of course, all of this became unnecessary when I bought my HP-35.
Michael

Re: Computing Square Roots
Message #6 Posted by Garth Wilson on 29 Mar 2009, 1:05 a.m.,
in response to message #5 by Michael de Estrada

We spent almost no time at all on the paired-digit method, so I doubt that any of us remembered it more than a few days. What you can do on paper though, in a pinch, is to estimate the answer, divide the original number by it, and then know that the real answer is between your estimate and the preliminary answer; so you can take the average of the two, skew it down slightly if you want to, and repeat. It's not as efficient, but it's unlikely you'll forget it. Back then, myself, I just used the slide rule.

Re: Computing Square Roots
Message #7 Posted by Marcus von Cube, Germany on 29 Mar 2009, 5:03 a.m.,
in response to message #6 by Garth Wilson

Your method reminds me on the starting point for Newtons method:

f(x)=x²-A f'(x)=2x x_n+1=x_n-f(x_n)/f'(x_n) x_n+1=x_n-(x_n²-A)/2x_n
Let's try an example: A=5
x₀=2 x₁=2-(2²-5)/(2*2)=2+1/4=2.25 x₂=2.25-(2.25²-5)/(2*2.25)=2.2361


Re: Computing Square Roots
Message #8 Posted by Khanh-Dang Nguyen Thu-Lam on 29 Mar 2009, 6:43 a.m.,
in response to message #7 by Marcus von Cube, Germany

Interestingly enough, the Newton's method is the same than the formula one gets with the geometrical interpretation. Let's rewrite the Newton's formula as:
x_n+1 = (x_n + A/x_n) / 2 .
At step n, x_n may be interpreted as the width of a rectangle of height A/x_n. The rectangle's area is just width times height, that is to say, A.
At step n+1, x_n+1 is the mean of the previous width and previous height. x_n+1 is the new width and A/x_n+1 is the new height. You still have a rectangle of area A, but it looks more and more like a square: the width is getting smaller, whereas the height is getting larger.
When n goes to infinity, the width and height converge to the same value which may be interpreted as the side of a square of area A.

Re: Computing Square Roots
Message #9 Posted by Thomas Klemm on 29 Mar 2009, 11:15 a.m.,
in response to message #7 by Marcus von Cube, Germany

The Horner schema can be used with Newton's method as well:
1 0 -5
2 1 2 -1 = f(2) 1 4 = f'(2) 1 -------------------------------------------------------------------- 1 4 -1 dx = 0.25 = - (-1) / 4
0.25 1 4.25 0.0625 1 4.5 1 -------------------------------------------------------------------- 1 4.5 0.0625 dx = - 0.01389 = - 0.0625 / 4.5
-0.013889 1 4.486111 0.000193 1 4.472222 1 --------- 2.236111 ========

Or you might use Cochran's trick:
5 0 -25
2 5 10 -5 5 20 5 -------------------------------------------------------------------- 5 20 -5 dx = 0.25 = - (-5) / 20
0.25 5 21.25 0.3125 5 22.5 5 -------------------------------------------------------------------- 5 22.5 0.3125 dx = - 0.01389 = - 0.3125 / 22.5
-0.013889 5 22.430556 0.000965 5 22.361111 5 =========
Now you get the sum of the differences gratis.

However I assume Newton would have used a different method as I already noted in a previous thread.
As can be seen from the recursion formula it's a cubic polynomial we have to evaluate but now the difference contains only a cheap division by 2:
1 - d * x_i² x_i+1 = x_i + x_i ----------- 2

Here's the same calculation using Horner schema:
-5 0 1 0
0.4 -5 -2 0.2 0.08 -5 -4 -1.4 -5 -6 -5 -------------------------------------------------------------------- 0.4 -5 -6 -1.4 0.08 dx = 0.04 = 0.08 / 2
0.04 -5 -6.2 -1.648 0.01408 -5 -6.4 -1.904 -5 -6.6 -5 -------------------------------------------------------------------- 0.44 -5 -6.6 -1.90400 0.01408 dx = 0.00704 = 0.01408 / 2
0.00704 -5 -6.6352 -1.95071 0.00035 -5 -6.6704 -1.99767 -5 -6.7056 -5 -------------------------------------------------------------------- 0.44704 -5 -6.7056 -1.99767 0.00035 dx = 0.00017 = 0.00035 / 2
0.00017 -5 -6.70647 -1.99883 0.00000 -5 -6.70733 -2.00000 -5 -6.70820 -5 -------------------------------------------------------------------- 0.44721 -5 -6.70820 -2.00000 0.00000
0.44721 * 5 ------- 2.23607 =======
Of course later in his life Newton used his HP-10C a time voyager left on his journey.

Edited: 29 Mar 2009, 11:24 a.m.

Re: Computing Square Roots
Message #10 Posted by John B. Smitherman on 29 Mar 2009, 3:52 p.m.,
in response to message #1 by Thomas Klemm

Thanks for the post Thomas. With the introduction of pocket calculators I hadn't thought about this for quite awhile but I remember learning the Babylonian method as discussed here:
http://en.wikipedia.org/wiki/Methods_of_computing_square_roots
Regards,
John

Re: Computing Square Roots
Message #11 Posted by cyrille de Br�bisson on 30 Mar 2009, 12:00 p.m.,
in response to message #1 by Thomas Klemm

hello,
the method used by all HO calculators from the 71 is described at: in the HP Solve: Volume 5 (June 2008) newsletter at: http://h20331.www2.hp.com/Hpsub/cache/580500-0-0-225-121.html
regards, cyrille

Re: Computing Square Roots
Message #12 Posted by Thomas Klemm on 30 Mar 2009, 3:26 p.m.,
in response to message #11 by cyrille de Br�bisson

These are the figures copied from the newsletter where the square root of 2 is calculated:

i mantissa result increment
12 100,000,000,000 100,000,000,000 100,000,000,000 11 200,000,000,000 140,000,000,000 10,000,000,000 10 595,000,000,000 141,000,000,000 1,000,000,000 9 302,000,000,000 141,400,000,000 100,000,000 8 191,800,000,000 141,420,000,000 10,000,000 7 503,795,000,000 141,421,000,000 1,000,000 6 795,315,500,000 141,421,300,000 100,000 5 882,088,750,000 141,421,350,000 10,000 4 335,606,320,000 141,421,356,000 1,000 3 527,636,078,000 141,421,356,200 100 2 103,372,009,355 141,421,356,230 10 1 437,705,999,155 141,421,356,237 1
And that's what you get when using the Horner schema:

5 0 -10
5 10 -5 5 20 10 5 100 -500
5 110 -395 5 120 -280 5 130 -155 5 140 -20 5 150 125 5 1400 -2000
5 1410 -595 5 1420 820 5 14100 -59500
5 14110 -45395 5 14120 -31280 5 14130 -17155 5 14140 -3020 5 14150 11125 5 141400 -302000
5 141410 -160595 5 141420 -19180 5 141430 122245 5 1414200 -1918000
5 1414210 -503795 5 1414220 910420 5 14142100 -50379500
5 14142110 -36237395 5 14142120 -22095280 5 14142130 -7953155 5 14142140 6188980 5 141421300 -795315500
5 141421310 -653894195 5 141421320 -512472880 5 141421330 -371051555 5 141421340 -229630220 5 141421350 -88208875 5 141421360 53212480 5 1414213500 -8820887500
5 1414213510 -7406673995 5 1414213520 -5992460480 5 1414213530 -4578246955 5 1414213540 -3164033420 5 1414213550 -1749819875 5 1414213560 -335606320 5 1414213570 1078607245 5 14142135600 -33560632000
5 14142135610 -19418496395 5 14142135620 -5276360780 5 14142135630 8865774845 5 141421356200 -527636078000
5 141421356210 -386214721795 5 141421356220 -244793365580 5 141421356230 -103372009355 5 141421356240 38049346880 5 1414213562300 -10337200935500
5 1414213562310 -8922987373195 5 1414213562320 -7508773810880 5 1414213562330 -6094560248555 5 1414213562340 -4680346686220 5 1414213562350 -3266133123875 5 1414213562360 -1851919561520 5 1414213562370 -437705999155 5 1414213562380 976507563220 5 14142135623700 -43770599915500
As I already mentioned: looking at the same thing from a different angle.
PS: Our first lines disagree. Could it be that the first line shold read instead:

12 500,000,000,000 100,000,000,000 100,000,000,000
At which line of the program should I assume is the printing of the variables?
If I suppose it's before the while-loop then result is still 0 wheras when it's after the while-loop mantissa has already been diminished.

[ Return to Index | Top of Index ]

Go back to the main exhibit hall