Login

Thomas Klemm · 03-27-2009, 05:43 PM

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

Bob Patton · 03-27-2009, 09:00 PM

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

Thomas Klemm · 03-28-2009, 07:24 AM

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

Karl Schneider · 03-28-2009, 08:43 PM

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

Michael de Estrada · 03-28-2009, 09:33 PM

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

Garth Wilson · 03-29-2009, 01:05 AM

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.

Marcus von Cube, Germany · 03-29-2009, 05:03 AM

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

Khanh-Dang Nguyen Thu-Lam · 03-29-2009, 06:43 AM

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.

Thomas Klemm · 03-29-2009, 11:15 AM

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.

John B. Smitherman · 03-29-2009, 03:52 PM

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

cyrille de Brébisson · 03-30-2009, 12:00 PM

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

Thomas Klemm · 03-30-2009, 03:26 PM

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.

Login
Username:
Password:	Lost Password?
	Remember me

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