Re: Program question Message #11 Posted by Thomas Klemm on 31 July 2007, 3:44 p.m., in response to message #1 by Vincze
Instead of calculating the square root of a number d the inverse of the square root may be calculated followed by the multiplication of d:
1
sqrt(d) = ------- * d
sqrt(d)
Use an initial value:
x0 ~ 1 / sqrt(d)
and calculate x1, x2, x3, ...
with the following formula:
1 - d * xi2
xi+1 = xi + xi -----------
2
This method has the advantage that only a division by 2 is needed.
Here's an HP-11c program:
01 LBL A 14 *
02 STO I 15 -
03 1 16 *
04 + 17 2
05 1/x 18 /
06 ENTER 19 RND
07 LBL 0 20 x#0
08 + 21 GTO 0
09 ENTER 22 +
10 x^2 23 RCL I
11 1 24 *
12 x<>y 25 RTN
13 RCL I
I guess it's easy to convert that into a program for the HP-35s.
Now let's suppose we want to calculate sqrt(5) only with paper
and pencil using an initial guess x0 = 0.4.
At the beginning the multiplications are easy to perform since we have only few significant digits whereas at the end the numbers become smaller and smaller.
However there's a slight difference compared to the program above.
Since xi+1 = xi + dx the square of the next value is calculated using the square of the previous:
(xi+1)2 = xi2 + (2*xi + dx)*dx
1.00000
.....................................................................
0.00000 ---> 0.00000 5.00000
+ 0.40000 ----------------------> * 0.40000
------- -------
0.40000 ---> + 0.40000 2.00000
======= -------
0.40000 ---> * 0.40000
-------
0.80000 ---> - 0.80000
-------
* 0.20000 <--------------------------------------------- 0.20000
-------
0.08000
/ 2
-------
0.04000
=======
.....................................................................
0.40000 ---> 0.40000 5.00000
+ 0.04000 ----------------------> * 0.04000
------- -------
0.44000 ---> + 0.44000 0.20000
======= -------
0.84000 ---> * 0.84000
-------
0.16800 ---> - 0.16800
-------
* 0.03200 <--------------------------------------------- 0.03200
-------
0.01408
/ 2
-------
0.00704
=======
.....................................................................
0.44000 ---> 0.44000 5.00000
+ 0.00704 ----------------------> * 0.00704
------- -------
0.44704 ---> + 0.44704 0.03520
======= -------
0.88704 ---> * 0.88704
-------
0.03122 ---> - 0.03122
-------
* 0.00078 <--------------------------------------------- 0.00078
-------
0.00035
/ 2
-------
0.00017
=======
.....................................................................
0.44704 ---> 0.44704 5.00000
+ 0.00017 ----------------------> * 0.00017
------- -------
0.44721 ---> + 0.44721 0.00087
======= -------
0.89425 ---> * 0.89425
-------
0.00078 ---> - 0.00078
-------
* 0.00000 <--------------------------------------------- 0.00000
-------
0.00000
/ 2
-------
0.00000
=======
.....................................................................
0.44721
* 5.00000
-------
2.23607
=======
Credits to Newton for this marvelous algorithm.
Edited: 31 July 2007, 4:09 p.m.
|