Threaded Mode | Linear Mode

Gerson W. Barbosa · 06-04-2021, 10:47 PM

Thank you very much, C.Ret, for the HP-15C version!
The HP-15C was my first calculator and I am supposed to know how to program it, yet I had no idea on how to start.
The results almost exactly match my latest table, except for 1, 2 or 3 ULPs. There was only one result with a 5-ULP difference, but it turned out to be still close to the actual result, but on the other side.
Only about 8 seconds on my HP-15C, purchased brand-new in ‘83 or ‘84. As a comparison, one AGM iteration takes 2.5 seconds, but then again my program might not be nearly as optimized as it could.

C.Ret · (This post was last modified: 06-05-2021 05:03 AM by C.Ret.)

Code:

001- ►LBL E

002-     STO 2  x²  x<>y  STO 1  x²  +  STO 3  1  STO 0        ; Initiate  R2:b  R1:a  R3:S◄a²+b²  R0:t

011-     ►LBL 0

012-        RCL 1  ENTER^  RCL*2  √x  x<> 2  -  2  STO*0  /    ; Compute   R2:b◄√ab  R0:t◄t+t     k◄(a-b)/2 

021-        STO-1  x²  RCL*0  STO-3                            ;           R1:a◄a-k            R3:S◄S-t*k²

025-        RCL+3  RCL 3  x<y?  GTO 0                          ; Repeat while  S+t*k² > S                (Test 8) 

029-    RCL/2  PI  *

032- RTN                                                       ; Return    P(a,b)◄S/B*π

ijabbott · 06-05-2021, 10:31 AM

I found this paper that uses a modified AGM (MAGM) sequence. It may be of interest to some people (expecially Albert Chan if he hasn't already seen it!).

An Eloquent Formula for the Perimeter of an Ellipse
Semjon Adlaj

C.Ret · (This post was last modified: 06-05-2021 03:41 PM by C.Ret.)

Thank a lot for the publication reference.

I was wandering what AGM stand for and I was at the point asking for explanations. Through this instructive reading, I have now the good definitions and explanations of what are AGM and MAGM. That all I need for a better understanding on how my program works !

Gerson W. Barbosa · 06-05-2021, 08:00 PM

(06-04-2021 06:12 PM)Gerson W. Barbosa Wrote: ———————

This will return exact results when y = 1:

p(a, b) ~ π(a - b)(y + 1/(4y - 1/(4y - 3/(4y - 3/(4y - 11/(12y - (4/(2y - 1) + (2592-825π)/(588-187π))))))))

Code:

… 25 C=4/(2*Y-1)+(2592-825*PI)/(588-187*PI) …

While this forces an exact result for y = 0, it only spoils the other results.
The same happens when the 1/3 constant is used.

A better option for both 10 and 12-digit results is to discard even the 1/3 constant:

p(a, b) ~ π(a - b)(y + 1/(4y - 1/(4y - 3/(4y - 3/(4y - 11/(12y - (4/(2y - 1))))))))

Code:

…

25 C=4/(2*Y-1)

…

In C.Ret’s program just delete the steps 16 through 18 [ 3 1/x + ].

Albert Chan · 06-05-2021, 10:35 PM

(06-05-2021 08:00 PM)Gerson W. Barbosa Wrote: A better option for both 10 and 12-digit results is to discard even the 1/3 constant:

p(a, b) ~ π(a - b)(y + 1/(4y - 1/(4y - 3/(4y - 3/(4y - 11/(12y - (4/(2y - 1))))))))

Try replacing last constant 1 as (19/18)/y.
This keep rel error below 1 ppm, well until h = ((a-b)/(a+b))^2 > 0.9616

C.Ret · (This post was last modified: 06-07-2021 07:25 PM by C.Ret.)

(06-05-2021 08:00 PM)Gerson W. Barbosa Wrote: A better option for both 10 and 12-digit results is to discard even the 1/3 constant:
p(a, b) ~ π(a - b)(y + 1/(4y - 1/(4y - 3/(4y - 3/(4y - 11/(12y - (4/(2y - 1))))))))
[...]
In C.Ret’s program just delete the steps 16 through 18 [ 3 1/x + ].

(06-05-2021 10:35 PM)Albert Chan Wrote: Try replacing last constant 1 as (19/18)/y. This keep rel error below 1 ppm, well until h = ((a-b)/(a+b))^2 > 0.9616

Code:

001- ►LBL E

002-    x=y ?  GTO 8                                ; Test for circle

004-    STO 0  x<>y  STO+0  -  CHS  STO 1  STO/0    ; Initiate    R0:y=(a+b)/(a-b)    R1:(a-b)  

011-    19  RCL/0  18  /  2  x²  LSTx  GSB 0        ; Initiate    c_ 0 _ = 4/(2*y - 19/18y )

021-    11  GSB.2  GSB 3  GSB 3  GSB 1  GSB 1       ; Iterate     c_i+1_ = n/(d*y -  c_i_  )

028-    RCL+0  RCL*1  GTO 9                         ; Preompute   (a-b)*(y+c)  

031-     ►LBL.2    12  GTO 0                        ; ── ─── ─── d=12 ─┐    

035-     ►LBL 1     1  GTO 4                        ; ── n=1 ─┐        │

038-     ►LBL 3     3                               ; ── n=3 ─┤        │

040-     ►LBL 4     4                               ;         └─ d=4 ──┤

042-     ►LBL 0        RCL 0  *  R^  -  /  RTN      ;                  └─ c'=n/(d*y-c) ───

049- ►LBL 8     2  *

052- ►LBL 9     PI  *                               ; Finalize   P(a,a)=2πa 

055- RTN                                            ;       or   P(a,b)~π(a-b)(y+1/(4y-1/(4y-3/(4y-3/(4y-11/(12y-(4/(2y-19/18y))))))))

                                                         ;                                                         with y=(a+b)/(a-b)

\( P(a,b)\approx\pi(a-b)(y+1/(4y-1/(4y-3/(4y-3/(4y-11/(12y-(4/(2y-19/18y)))))))) \)

\( P(a,b)\approx\pi(a-b)\frac{110592y^8-122880y^6+18704y^4+7488y^2-627}{110592y^7-150528y^5+54608y^3-4244y} \) with \( y=\frac{a+b}{a-b} \)

\( P(a,b)\approx\pi(a+b)\frac{110592-122880h+18704h^2+7488h^3-627h^4}{110592-150528h+54608h^2-4244h^3} \) with \( h=\left(\frac{a-b}{a+b}\right)^2 \)

Code:

a‌‌‌‌‌‏   b     P(a,b)               Calc'Lap

20   20      125,6637062      1"54

20   19      122,5422527      9"03

20   10      96,88448221      8"83

20    5      85,78421626      9"48

20    1      80,38844081      9"11

20    0,05   80,00067833      9"04

20    0,01   79,99830115      8"99

20    0      79,99793952      8"79

7 DATA 19,18, 4,2, 11,12, 3,4, 3,4, 1,4, 1,4
10 INPUT "Elps Radii A,B ";A,B @ IF A<>B THEN Y=(A+B)/(A-B) ELSE P=2*PI*A @ GOTO 30
20 C=0 @ FOR I=1 TO 7 @ READ N,D @ C=N/(D*Y-C) @ NEXT I @ P=PI*(A-B)*(Y+C)
30 DISP USING 40;A,B,P
40 IMAGE "P("K",",K")=",4D.8D

Gerson W. Barbosa · (This post was last modified: 07-14-2021 11:36 PM by Gerson W. Barbosa.)

(06-05-2021 12:31 AM)C.Ret Wrote:

Code:

001- ►LBL E 002- STO 2 x² x<>y STO 1 x² + STO 3 1 STO 0 ; Initiate R2:b R1:a R3:S◄a²+b² R0:t 011- ►LBL 0 012- RCL 1 ENTER^ RCL*2 √x x<> 2 - 2 STO*0 / ; Compute R2:b◄√ab R0:t◄t+t k◄(a-b)/2 021- STO-1 x² RCL*0 STO-3 ; R1:a◄a-k R3:S◄S-t*k² 025- RCL+3 RCL 3 x<y? GTO 0 ; Repeat while S+t*k² > S (Test 8) 029- RCL/2 PI * 032- RTN ; Return P(a,b)◄S/B*π

Just for the record, here is a version of your program using three registers, by former forum member Mike (Stgt), from Jun/06. Actually, he wrote it for the HP-15C, like you did, but I prefer the HP-32S II, because it’s faster.

Code:

P01 LBL P

P02 STO B

P03 x²

P04 x<>y

P05 STO A

P06 x²

P07 +

P08 1

P08 STO C

Q01 LBL Q

Q02 RCL A

Q03 ENTER

Q04 RCL× B

Q05 SQRT

Q06 x<> B

Q07 -

Q08 2

Q09 STO× C

Q10 ÷

Q11 STO— A

Q12 x²

Q13 RCL× C

Q14 x<>y

Q15 R↓

Q16 -

Q17 x<>y

Q18 x>y?

Q19 GTO Q

Q20 RCL÷ B

Q21 π

Q22 ×

Q23 RTN

P: CK=F0E5 013.5 bytes

Q: CK=D8B7 034.5 bytes

Regards,

Gerson.

floppy · 04-18-2023, 09:24 AM

(05-26-2021 11:49 AM)MeindertKuipers Wrote: Some background on the approximations for the perimeter of an ellipse, amazing video to watch (like most of this series)

Remark: Parker is an "non-reference" (a bit a tik-tok math-touch video. however, took me time to come to this conclusion).
Comments regarding this video page 8 of https://indico-hlit.jinr.ru/event/187/co...SAdlaj.pdf