Computing Square Roots



#13

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:

x2 - 41 = 0.

First we apply a trick which can be found in Cochran's patent and multiply the whole equation by 5:

5x2 + 0x - 205 = 0

Now let's consider the following simple coordinate transformations:

  1. shift: x' = x - 1
  2. zoom: x' = 10 * x

What happens to a quadratic equation if we apply these transformations?

Shift

P(x) = ax2 + bx + c
= Q1(x)(x-1) + c'
= (Q2(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

ax2 + bx + c = 0    | * 100
100ax2 + 100bx + 100c = 0
ax'2 + 10bx' + 100c = 0
a'x'2 + 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 x2 = 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


#14

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


#15

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

#16

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


#17

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


#18

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.


#19

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

f(x)=x2-A
f'(x)=2x
xn+1=xn-f(xn)/f'(xn)
xn+1=xn-(xn2-A)/2xn
Let's try an example: A=5
x0=2
x1=2-(22-5)/(2*2)=2+1/4=2.25
x2=2.25-(2.252-5)/(2*2.25)=2.2361


#20

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:

xn+1 = (xn + A/xn) / 2 .

At step n, xn may be interpreted as the width of a rectangle of height A/xn. The rectangle's area is just width times height, that is to say, A.

At step n+1, xn+1 is the mean of the previous width and previous height. xn+1 is the new width and A/xn+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.

#21

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 * xi2
xi+1 = xi + xi -----------
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.

#22

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

#23

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


#24

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.


Possibly Related Threads...
Thread Author Replies Views Last Post
  [HP-Prime CAS] "Warning, ^ (Command) Is ambiguous on non square matrices"?? CompSystems 1 366 12-07-2013, 07:15 PM
Last Post: CompSystems
  HP prime: logs and roots Alberto Candel 5 384 11-20-2013, 12:31 PM
Last Post: Manolo Sobrino
  Computing pi with the PC-1300S Kiyoshi Akima 0 190 11-17-2013, 12:24 AM
Last Post: Kiyoshi Akima
  A brand new calculator benchmark: "middle square method seed test" Pier Aiello 25 1,417 09-13-2013, 01:58 PM
Last Post: Pier Aiello
  Square Root Simplifier for HP39gII Mic 4 368 03-11-2013, 08:25 AM
Last Post: Eddie W. Shore
  Books on computing Elementary Functions Nick_S 4 270 09-14-2012, 07:33 AM
Last Post: David Hayden
  [WP34S] Improving the Upper-Tail F and Chi-Square Les Wright 7 452 05-12-2012, 06:31 PM
Last Post: Paul Dale
  [WP34S] Chi-Square Functions Broken in 2935 Les Wright 7 493 05-03-2012, 11:10 AM
Last Post: Les Wright
  HP-51G: Back to the honorable roots! Martin Paech 24 1,252 04-08-2012, 01:39 AM
Last Post: db (martinez, ca.)
  35s - find roots of 3rd and higher order equations Chris C 11 606 02-25-2012, 02:55 PM
Last Post: Thomas Klemm

Forum Jump: