Post Reply 
Perimeter of the Ellipse (HP-15C)
06-04-2021, 10:47 PM
Post: #21
RE: Perimeter of the Ellipse (HP-15C)
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.
Find all posts by this user
Quote this message in a reply
06-05-2021, 12:31 AM (This post was last modified: 06-05-2021 05:03 AM by C.Ret.)
Post: #22
RE: Perimeter of the Ellipse (HP-15C)
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*π
Find all posts by this user
Quote this message in a reply
06-05-2021, 10:31 AM
Post: #23
RE: Perimeter of the Ellipse (HP-15C)
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

— Ian Abbott
Find all posts by this user
Quote this message in a reply
06-05-2021, 03:40 PM (This post was last modified: 06-05-2021 03:41 PM by C.Ret.)
Post: #24
RE: Perimeter of the Ellipse (HP-15C)
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 !
Find all posts by this user
Quote this message in a reply
06-05-2021, 08:00 PM
Post: #25
RE: Perimeter of the Ellipse (HP-15C)
(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 + ].
Find all posts by this user
Quote this message in a reply
06-05-2021, 10:35 PM
Post: #26
RE: Perimeter of the Ellipse (HP-15C)
(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
Find all posts by this user
Quote this message in a reply
06-06-2021, 07:07 AM (This post was last modified: 06-07-2021 07:25 PM by C.Ret.)
Post: #27
RE: Perimeter of the Ellipse (HP-15C)
(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
Find all posts by this user
Quote this message in a reply
07-14-2021, 11:34 PM (This post was last modified: 07-14-2021 11:36 PM by Gerson W. Barbosa.)
Post: #28
RE: Perimeter of the Ellipse (HP-15C)
(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.
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)