Perimeter of the Ellipse (HP-15C)
06-04-2021, 10:47 PM
Post: #21
 Gerson W. Barbosa Senior Member Posts: 1,404 Joined: Dec 2013
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.
06-05-2021, 12:31 AM (This post was last modified: 06-05-2021 05:03 AM by C.Ret.)
Post: #22
 C.Ret Member Posts: 122 Joined: Dec 2013
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*π
06-05-2021, 10:31 AM
Post: #23
 ijabbott Senior Member Posts: 1,023 Joined: Jul 2015
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

— Ian Abbott
06-05-2021, 03:40 PM (This post was last modified: 06-05-2021 03:41 PM by C.Ret.)
Post: #24
 C.Ret Member Posts: 122 Joined: Dec 2013
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 !
06-05-2021, 08:00 PM
Post: #25
 Gerson W. Barbosa Senior Member Posts: 1,404 Joined: Dec 2013
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 + ].
06-05-2021, 10:35 PM
Post: #26
 Albert Chan Senior Member Posts: 1,591 Joined: Jul 2018
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
06-06-2021, 07:07 AM (This post was last modified: 06-07-2021 07:25 PM by C.Ret.)
Post: #27
 C.Ret Member Posts: 122 Joined: Dec 2013
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
07-14-2021, 11:34 PM (This post was last modified: 07-14-2021 11:36 PM by Gerson W. Barbosa.)
Post: #28
 Gerson W. Barbosa Senior Member Posts: 1,404 Joined: Dec 2013
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.
 « Next Oldest | Next Newest »

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