(HP41) Iterated function with convergence to square root (Heron method)

[floppy]

Since RPN and FOCAL with it stack is efficient for iterated functions, here a posting of the iterated function of aN=(an + b/an) / 2 which result will be SQRT(b).
Other way to represent this is f(x,y)=(x+y/x,y). Then f(f(f(..(a0,b))))=(sqrt(b),b).
The start value a0 can be anything (result will have the negativ sign if the start value a0<0).
Nothing really new:
- no register used in comparison to the programm available in page 11 in https://literature.hpcalc.org/community/...ram-de.pdf )
- start value can be defined (will be defined as b if a0=0)
- no calculation done if b=0

Code:

; SQRT as iterated function, "möbius transform like" = only + * or /  
;
; Execution/Inputs
;  a0 ENTER b XEQ ALPHA SQRT ALPHA or
;  b XEQ ALPHA SQRT ALPHA
;  ... b number which should be squared
;  ... a0 is a start value for the iteration
;  ... if a0 is negativ, the result will be negativ value of the square root
;
; Outputs
;  results in stack X: sqrt(b)
;
; Modules used
;  None
;
; use Register none
;
; stack effect: all overwritten (X Y Z LASTX)
;
; create raw file with "rpncomp --raw-output SQRT.TXT", upload
;  in V41 emulator and store into virtual drive with "ilper"
;  .. hp41uc https://sourceforge.net/projects/hp41uc/ or
;     rpncomp https://github.com/hth313/Calypsi-tool-chains
;  .. V41   https://hp.giesselink.com/v41.htm
;  .. ilper https://hp.giesselink.com/hpil.htm
; 
; under CC BY SA CreativeCommons 4.0 
;  floppy @ https://www.hpmuseum.org/forum/
;
; idea taken from here
; https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Heron's_method
;
; change log
;  date 2023 04 24 creation
;  date 2023 04 25 testing..
;
; comments
;  aN=(an + b/an) / 2
;  .. a number to be squared.. it mean aN (N-> infinite) = sqrt(a)
;
LBL "SQRT"
;                  X            Y           Z          T
;                  b           a0
X<=0?        ; in case b <= 0 dont do anything
RTN          ; ... just go out 
X<>Y         ;     a0           b  
X=0?         ; in case a = 0 set a to b
XEQ 00
X<>Y         ;     b            a0
;
STO Z        ;     b            a0          b
X<>Y         ;     a0           b          b
LBL 01
  RCL Z      ;     b            an          XX          b
  X<>Y       ;     an           b           XX          b
  STO Z      ;     an           b           an          b
  /          ;     b/an         an          b           b
  LASTX      ;     an          b/an         an          b    
  +          ;  (an+b/an)       an          b           b
  2
  /          ; (an+b/an)/2=aN   an          b           b
  X#Y?
  GOTO 01
RTN
LBL 00       ;     X            Y           Z          T
  RDN        ;     b            
  ENTER      ;     b            a(=b)         0     
RTN      
END

[C.Ret]

Hi,

It reminds me of a program, using this exact method, written for a silly challenge on a French forum. (Look for MPO n°46 on Silicium).


001 LBL"HERON"  RCL Y  RCL Y  /  +  2  /  END  

Usage:
Enter the number b whose square root you want to determine.
Press ENTER^.
Optionally enter the approximation a0.
Run the XEQ (ALPHA)HERON(ALPHA) program.

A better approximation is displayed.
Repeat the computation by pressing R/S . You will see the successive approximations quickly converging towards the square root of b.

To extract another root, simply enter a new value, press ENTER^, enter an optional guess approximation and restart with R/S.


Example:
2    FIX 9                     2,000000000
1.5  XEQ (ALPHA)HERON(ALPHA)   1,416666667
     R/S                       1,414215687
     R/S                       1,414213563
     R/S                       1,414213563
     ...                                  

[Thomas Klemm]

(04-25-2023 03:21 PM)floppy Wrote:  - no register used in comparison to the program available in page 11 in https://literature.hpcalc.org/community/...ram-de.pdf )

Why not compare it to the program on the next page, which also doesn't use registers?

[floppy]

One more https://www.youtube.com/watch?v=uMCWgcBv...ndreDumont

[floppy]

(04-25-2023 06:29 PM)Thomas Klemm Wrote:  

(04-25-2023 03:21 PM)floppy Wrote:  - no register used in comparison to the program available in page 11 in https://literature.hpcalc.org/community/...ram-de.pdf )

Why not compare it to the program on the next page, which also doesn't use registers?

correct. I missed this page 12. However, the program is with 2 variables = a bit different (the start value is something I will consider for another topic).

[Thomas Klemm]

(04-25-2023 05:51 PM)C.Ret Wrote:  written for a silly challenge

A bit shorter, albeit sillier:

Initialisation

CLΣ
ENTER
ENTER
ENTER

Iteration

/
Σ+
MEAN

[Thomas Klemm]

(04-25-2023 06:30 PM)floppy Wrote:  One more (…)

HP-41C:

Code:

STO 00
SIGN
LBL 00
RCL 00
RCL Y
/
LASTX
+
2
/
X#Y?
GTO 00
END

HP-42S:

Code:

00 { 16-Byte Prgm }
01 STO 00
02 SIGN
03▸LBL 00
04 RCL 00
05 RCL÷ ST Y
06 RCL+ ST Y
07 2
08 ÷
09 X≠Y?
10 GTO 00
11 END