(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