Perimeter of the Ellipse (HP-15C)

[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