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

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

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
...

(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?

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

(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).

(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

(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