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:
x_{0} ~ 1 / sqrt(d)
and calculate x_{1}, x_{2}, x_{3}, ...
with the following formula:
1  d * x_{i}^{2}
x_{i+1} = x_{i} + x_{i} 
2
This method has the advantage that only a division by 2 is needed.
Here's an HP11c 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 HP35s.
Now let's suppose we want to calculate sqrt(5) only with paper
and pencil using an initial guess x_{0} = 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 x_{i+1} = x_{i} + dx the square of the next value is calculated using the square of the previous:
(x_{i+1})^{2} = x_{i}^{2} + (2*x_{i} + 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.
