Perimeter of the Ellipse (HP-15C)

[Gerson W. Barbosa]

Code:


   001 {    42 21 15 } f LBL E
   002 {       44  1 } STO 1
   003 {           2 } 2
   004 {          20 } ×
   005 {          34 } x↔y
   006 {       44  2 } STO 2
   007 {    45 40  1 } RCL + 1
   008 {    43 30  6 } g TEST x≠y
   009 {          10 } ÷
   010 {           1 } 1
   011 {          30 } -
   012 {       43 11 } g x²
   013 {           3 } 3
   014 {           2 } 2
   015 {          34 } x↔y
   016 {          10 } ÷
   017 {       43 36 } g LSTx
   018 {          48 } .
   019 {           3 } 3
   020 {           6 } 6
   021 {           1 } 1
   022 {       43  2 } g →H
   023 {          20 } ×
   024 {          30 } -
   025 {           4 } 4
   026 {          30 } -
   027 {          15 } 1/x
   028 {           1 } 1
   029 {          48 } .
   030 {           5 } 5
   031 {          40 } +
   032 {    42  4  1 } f Χ↔ 1
   033 {       45  1 } RCL 1
   034 {          14 } y^x
   035 {       45  2 } RCL 2
   036 {       43 36 } g LSTx
   037 {          14 } y^x
   038 {          40 } +
   039 {           2 } 2
   040 {          10 } ÷
   041 {       45  1 } RCL 1
   042 {          15 } 1/x
   043 {          14 } y^x
   044 {           2 } 2
   045 {          20 } ×
   046 {       43 26 } g π
   047 {          20 } ×
   048 {       43 32 } g RTN

# ------------------------------------------------------------------------------
# Perimeter of the Ellipse
# ------------------------------------------------------------------------------
# Usage:
#
#   a ENTER b f E
# ------------------------------------------------------------------------------
# Formula:
#
#   p ≈ 2π[(aˢ + bˢ)/2]¹ᐟˢ
#
#   where
#
#   s = 3/2 + 1/(32/h² - 217h²/360 - 4)
#
#   and 
#
#   h = [(a - b)/(a + b)]
#
#  See https://www.hpmuseum.org/forum/post-126478.html#pid126478
#
# ------------------------------------------------------------------------------

# ------------------------------------------------------------------------------
# Listing generated by 
# HEWLETT·PACKARD 15C Simulator program
# Created with version 4.3.00
# ------------------------------------------------------------------------------
# © 2021 Torsten Manz
# http://hp-15c.homepage.t-online.de/download.htm
# ------------------------------------------------------------------------------


# --------

+-----+-----+-------+-------------+-------------+
|  a  |  b  |   h   |   result    |    exact    |
+-----+-----+-------+-------------+-------------+
| 20  | 20  | 0.000 | 125.6637062 | 125.6637061 |
| 20  | 19  | 0.026 | 122.5422527 | 122.5422527 |
| 20  | 18  | 0.053 | 119.4632087 | 119.4632087 |
| 20  | 17  | 0.081 | 116.4300496 | 116.4300496 |
| 20  | 16  | 0.111 | 113.4466716 | 113.4466716 | 
| 20  | 15  | 0.143 | 110.5174609 | 110.5174608 |
| 20  | 14  | 0.176 | 107.6473797 | 107.6473796 | 
| 20  | 13  | 0.212 | 104.8420720 | 104.8420720 |
| 20  | 12  | 0.250 | 102.1079954 | 102.1079955 | 
| 20  | 11  | 0.290 | 99.45258805 | 99.45258801 |
| 20  | 10  | 0.333 | 96.88448221 | 96.88448221 | 
| 20  |  9  | 0.379 | 94.41378489 | 94.41378488 |
| 20  |  8  | 0.429 | 92.05245051 | 92.05245038 | 
| 20  |  7  | 0.481 | 89.81479201 | 89.81479146 |
| 20  |  6  | 0.538 | 87.71820372 | 87.71820139 |
| 20  |  5  | 0.600 | 85.78422747 | 85.78421775 | 
| 20  |  4  | 0.667 | 84.04021830 | 84.04017816 |
| 20  |  3  | 0.739 | 82.52217428 | 82.52200588 | 
| 20  |  2  | 0.818 | 81.28023133 | 81.27948360 |
| 20  |  1  | 0.905 | 80.39234813 | 80.38851238 | 
| 20  |  0  | 1.000 | 80.03695848 | 80.00000000 | 
+-----+-----+-------|-------------+-------------+


Edited to fix a silly mistake in usage.

[OlidaBel]

Interesting subject , thanks for the post.
How did you get the exact value ? A computed integral (on your 15C off course) ?

[Gerson W. Barbosa]

I checked the results against WolframAlpha as I did eight years ago:

https://www.hpmuseum.org/cgi-sys/cgiwrap...ead=244138

Considering the integral evaluations took about 45 seconds on the significantly faster 42S, I won’t even try it on my old 15C (perhaps I’ll try it on the 15C LE later).

The approximation program takes almost seven seconds on my old15C.

Regards,

Gerson.

[MeindertKuipers]

Some background on the approximations for the perimeter of an ellipse, amazing video to watch (like most of this series)

[Albert Chan]

(01-19-2020 11:00 PM)Albert Chan Wrote:  \(\large p(a,b) = 2× \left( p({a+b \over 2}, \sqrt{ab}) - {\pi a b \over AGM(a,b)}\right)\)

In other words, calculate ellipse perimeter from a much less eccentric ellipse.

Code:

from math import sqrt, pi
def p(a, b, pi2=2*pi):
    'perimeter of ellipse, a>0, b>0, order does not matter'
    ab = a*b
    r = sqrt(ab)
    if r==b: return pi2*r, pi2/r
    p2, k = p((a+b)/2, r)
    return 2*p2-k*ab, k     # perimeter, 2*pi/agm(a,b)

>>> p(20,10)
(96.884482205476857, 0.43130312949992861)

With AGM, convergence is quadratic:

p(20, 10)
→ p(15.0, 14.142135623730951)
→ p(14.571067811865476, 14.564753151219703)
→ p(14.56791048154259, 14.567910139395549)
→ p(14.56791031046907, 14.567910310469069)

Last ellipse is practically a circle, thus return 2πr ≈ 91.532880019049259
Code then unwind back to the top, for p(20, 10) ≈ 96.884482205476857

Edit: code return (perimeter, 2*pi/agm(a,b))

>>> m = 0.5 # elliptic parameter m, 0 < m < 1
>>> [x/4 for x in p(1, sqrt(1-m))] # E(m), K(m)
[1.3506438810476749, 1.8540746773013719]

This may be why Matlab have ellipke that return both kinds of elliptic integrals.

[ttw]

The AGM is nice. Elliptic integrals (and their relatives) are a bit complicated.

[Gerson W. Barbosa]

(05-25-2021 10:05 PM)Gerson W. Barbosa Wrote:  

Code:


   001 {    42 21 15 } f LBL E
   002 {       44  1 } STO 1
   003 {           2 } 2
   004 {          20 } ×
   005 {          34 } x↔y
   006 {       44  2 } STO 2
   007 {    45 40  1 } RCL + 1
   008 {    43 30  6 } g TEST x≠y
   009 {          10 } ÷
   010 {           1 } 1
   011 {          30 } -
   012 {       43 11 } g x²
   013 {           3 } 3
   014 {           2 } 2
   015 {          34 } x↔y
   016 {          10 } ÷
   017 {       43 36 } g LSTx
   018 {          48 } .
   019 {           3 } 3
   020 {           6 } 6
   021 {           1 } 1
   022 {       43  2 } g →H
   023 {          20 } ×
   024 {          30 } -
   025 {           4 } 4
   026 {          30 } -
   027 {          15 } 1/x
   028 {           1 } 1
   029 {          48 } .
   030 {           5 } 5
   031 {          40 } +
   032 {    42  4  1 } f Χ↔ 1
   033 {       45  1 } RCL 1
   034 {          14 } y^x
   035 {       45  2 } RCL 2
   036 {       43 36 } g LSTx
   037 {          14 } y^x
   038 {          40 } +
   039 {           2 } 2
   040 {          10 } ÷
   041 {       45  1 } RCL 1
   042 {          15 } 1/x
   043 {          14 } y^x
   044 {           2 } 2
   045 {          20 } ×
   046 {       43 26 } g π
   047 {          20 } ×
   048 {       43 32 } g RTN

# ------------------------------------------------------------------------------
# Perimeter of the Ellipse
# ------------------------------------------------------------------------------
# Usage:
#
#   a ENTER b f E
# ------------------------------------------------------------------------------
# Formula:
#
#   p ≈ 2π[(aˢ + bˢ)/2]¹ᐟˢ
#
#   where
#
#   s = 3/2 + 1/(32/h² - 217h²/360 - 4)
#
#   and 
#
#   h = [(a - b)/(a + b)]
#
#  See https://www.hpmuseum.org/forum/post-126478.html#pid126478
#
# ------------------------------------------------------------------------------

# ------------------------------------------------------------------------------
# Listing generated by 
# HEWLETT·PACKARD 15C Simulator program
# Created with version 4.3.00
# ------------------------------------------------------------------------------
# © 2021 Torsten Manz
# http://hp-15c.homepage.t-online.de/download.htm
# ------------------------------------------------------------------------------


# --------

Please replace steps 3 through 15 with these:

Code:


   003 {          34 } x↔y
   004 {       44  2 } STO 2
   005 {          30 } -
   006 {       45  1 } RCL 1
   007 {    45 40  2 } RCL + 2
   008 {          10 } ÷
   009 {       43 11 } g x² 
   010 {           3 } 3
   011 {           2 } 2
   012 {          34 } x↔y
   013 {       43 20 } g x=0
   014 {       42  0 } g x!

Now 0 ENTER 0 f E will return “Error 0” (no problem, we all know the result is 0). On the other hand 0.5 ENTER 0.5 f E will return 3.141592654 comme il faut, instead of a division by zero error.

[PedroLeiva]

I have develop a program for HP67 to calculate perimeter, surface and volume of an ellipse. The inputs are diameters, but for calculation uses radius. When the 3 diameters are the same figure the body transforms into an sphere. The input and outputs can be printed or screen flashes (HP97/67), but also the results remains in the stack
In the attach PDF (3 pages A4) includes: a) program listing, b) formulas, c) Instructions, d) PGRM Crd diagram, e) example and calculation comparison for sphere using ellipse program (differences in the 8th. decimal only, dif.=+-1)

[Gerson W. Barbosa]

(05-27-2021 06:02 PM)PedroLeiva Wrote:  (differences in the 8th. decimal only, dif.=+-1)

Have you tested it for ellipses with h greater than 0.6? Although that statement should be valid for the extreme case in my examples (a = 20 and b = 0), the errors would show up starting from the fifth and sixth decimals (a =20; b = 19 and b = 18, respectively - anyway, that’s one more exact digit when compared to the approximation I have used, for these two examples).
You are using Cantrell-Ramanujan approximation, described at the end of this article.

Regards,

Gerson.

[PedroLeiva]

(05-27-2021 09:59 PM)Gerson W. Barbosa Wrote:  

(05-27-2021 06:02 PM)PedroLeiva Wrote:  (differences in the 8th. decimal only, dif.=+-1)

Have you tested it for ellipses with h greater than 0.6? Although that statement should be valid for the extreme case in my examples (a = 20 and b = 0), the errors would show up starting from the fifth and sixth decimals (a =20; b = 19 and b = 18, respectively - anyway, that’s one more exact digit when compared to the approximation I have used, for these two examples).
You are using Cantrell-Ramanujan approximation, described at the end of this article.

Regards,

Gerson.

I am using Ramanujan II-Cantrell: here a & b are radius (in my PDF the input is diameter)
H= [(a- b) / (a + b)]^2
P= π * (a+b) * [ 1 + 3H / (10+ SQRT(4-3H)) + (4/π - 14/11) * H^12]

For ellipses with H>0.6, here the examples:
a= 20 ----H= 0.546313800
b= 3 ----P= 82.52178335

a= 20 ---H= 0.669421488
b= 2 ---P= 81.27883093

a= 20 ---H= 0.904818560
b= 0.5 ---P= 80.11412754

For the other combinations of the radius, the results are:
a= 20 ---H= 1
b= 0 ---P= 80

a= 20 ---H= 0.0006557462
b= 19 ---P= 122.5422527

a= 20 ---H= 0.00277083
b= 18 ---P= 119.4632087

Please let me know your conclusions,
Pedro

[Gerson W. Barbosa]

(05-27-2021 11:28 PM)PedroLeiva Wrote:  I am using Ramanujan II-Cantrell: here a & b are radius (in my PDF the input is diameter)
H= [(a- b) / (a + b)]^2
P= π * (a+b) * [ 1 + 3H / (10+ SQRT(4-3H)) + (4/π - 14/11) * H^12]

For ellipses with H>0.6, here the examples:
a= 20 ----H= 0.546313800
b= 3 ----P= 82.52178335

a= 20 ---H= 0.669421488
b= 2 ---P= 81.27883093

a= 20 ---H= 0.904818560
b= 0.5 ---P= 80.11412754

For the other combinations of the radius, the results are:
a= 20 ---H= 1
b= 0 ---P= 80

a= 20 ---H= 0.0006557462
b= 19 ---P= 122.5422527

a= 20 ---H= 0.00277083
b= 18 ---P= 119.4632087

Please let me know your conclusions,

I was thinking of h = [(a - b) / (a + b)], not h = [(a - b) / (a + b)]^2. That is, I was actually interested in ellipses with h greater than 0.36, not 0.6.

Ramanujan-Cantrell is more accurate than Ramanujan-II, but only for very eccentric ellipsis (h close to 1).

I should have included a direct comparison between Ramanujan-Cantrell and Modified Muir approximations, but the spreadsheets below might give an idea on how these three approximations perform when compared to the exact results.

[PedroLeiva]

In my last reply something was not completely correct. Willy Kunz let me know that there was a mix of Ellipse and Ellipsoid concepts in calculation. For example, Perimeter and Area belongs to Ellipse; surface (not calculated before) and volume to Ellipsoid. So this new version of the program considers both: a) bottom row Crd for data and ellipse calculations, b) upper row Crd for ellipsoid calculations.

Pedro

[Gerson W. Barbosa]

Another method (under test)

p ~ π(a - b)(y + 1/(4y - 1/(4y - 3/(4y - 3/(4y - 4/(4y))))))

where

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

HP-75C, HP-71B program:

Code:


10 OPTION BASE 1
15 DIM N(5) 
20 INPUT A,B
25 IF A=B THEN B=-B @ Y=1 @ C=0 @ GOTO 65
30 Y=(A+B)/(A-B)
35 D=4*Y
40 C=0
45 FOR I=1 TO 5
50 READ N(I)
55 C=N(I)/(D-C)
60 NEXT I
65 P=PI*(A-B)*(Y+C)
70 DISP P
75 END
99 DATA 4,3,3,1,1

An HP-15C version could be slightly more accurate and perhaps faster than my first program.

[C.Ret]

Hello,

Great idea, but I am a bit lazy this morning to type so many lines of program on my quite full of stuff HP-71B.

\( P(a,b)=\pi(a-b)(y+\frac{1}{4y-\frac{1}{4y-\frac{3}{4y-\frac{3}{4y-\frac{4}{4y-\frac{3}{4y-\frac{2}{4y-\frac{2}{4y-\frac{5}{4y-0}}}}}}}}}) \) with \( y=\frac{a+b}{a-b} \)

Here is a version where coefficients don't have to be copied in an array; existing in program as DATA at once is enough and spare memory :

HP-71B version:

9 DATA 5,2,2,3,4,3,3,1,1
10 INPUT "ELIPSE 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 9 @ READ N @ C=N/(4*Y-C) @ NEXT I @ P=PI*(A-B)*(Y+C)
30 DISP USING 40;A,B,P
40 IMAGE "P("K",",K")=",4D.6D

Typical sample value:
P(20,20)= 125.663706
P(20,19)= 122.542253
P(19,20)= 122.542253
P(20,10)= 96.884482
P(20,1)= 80.388512

P(20,.2)= 80.016633

P(20,.1)= 79.998004
P(20,0)= 79.986832

HP-15C / HP-11C version:

Code:

001- ►LBL E
002-   STO 0  x↔y  STO+0  x=y?  GTO 1
007-   -  STO 1  STO/0  RCL-1
011-   5  GSB 0  2  GSB 0  2  GSB 0  3  GSB 0  4  GSB 0  3  GSB 0  3  GSB 0  1  GSB 0  1  GSB 0  RCL+0  RCL*1  GTO 2

032- ►LBL 0  4  RCL 0  *  R^  -  /  RTN

040- ►LBL 1  2  *
043- ►LBL 2  PI  *
046- RTN

Usage : a [ENTER^] b [ f ][ E ] returns ellipse perimeter P.

Registers:
R0: \( -y=\frac{a+b}{b-a} \)
R1: \( (b-a) \)

Edited several times to correct english, missing sentences and code typos.

[Gerson W. Barbosa]

Thank you, C.Ret, for the optimization.

It appears the actual approximation should be

π(a - b)(y + 1/(4y - 1/(4y - 3/(4y - 3/(4y - 11/(12y - 149/(57y)))))))

where y = (a + b)/(a - b)

Code:


10 DATA 11/(3-149/(228*Y*Y)),3,3,1,1
15 INPUT A,B
20 IF A<>B THEN Y=(A+B)/(A-B) ELSE P=2*PI*A @ GOTO 60
25 C=0
30 D=4*Y
35 FOR I=1 TO 5
40 READ N
45 C=N/(D-C)
50 NEXT I
55 P=PI*(A-B)*(Y+C)
60 DISP P

P(20, 0) = 79.9786371734

More terms would help or perhaps an optimization of the numerators so that the lower half of the table becomes more accurate.

P. S.: The error in the length of the orbit of Halley’s Comet is about only 117 km.

P(2667950000, 678281900) = 11 464 318 867.1 km
( actual result with these data: 11 464 318 984.1 km )

P. P. S.: Unlike the previous program, the following will work also on the HP-75C.

Code:


10 DATA 11,12,3,4,3,4,1,4,1,4
15 INPUT A,B
20 IF A<>B THEN Y=(A+B)/(A-B) ELSE P=2*PI*A @ GOTO 55
25 C=149/(57*Y)
30 FOR I=1 TO 5
35 READ N,D
40 C=N/(D*Y-C)
45 NEXT I
50 P=PI*(A-B)*(Y+C)
55 DISP P

[Albert Chan]

A simple formula that gives perimeter of ellipse with max rel error of 7ppm
h = ((a-b)/(a+b))^2

\(\displaystyle P(a,b) ≈ \pi (a+b) \left(
\frac{.001312}{1.2534-h} +
\frac{.85111}{3.8631-h} +
.192145\,h + .77864
\right) \)

[OlidaBel]

(06-03-2021 09:43 AM)C.Ret Wrote:  HP-15C / HP-11C version:

Code:

001- ►LBL E
002-   STO 0  x↔y  STO+0  x=y?  GTO 1
007-   -  STO 1  STO/0  RCL-1
011-   5  GSB 0  2  GSB 0  2  GSB 0  3  GSB 0  4  GSB 0  3  GSB 0  3  GSB 0  1  GSB 0  1  GSB 0  RCL+0  RCL*1  GTO 2

032- ►LBL 0  4  RCL 0  *  R^  -  /  RTN

040- ►LBL 1  2  *
043- ►LBL 2  PI  *
046- RTN

Usage : a [ENTER^] b [ f ][ E ] returns ellipse perimeter P.

Registers:
R0: \( -y=\frac{a+b}{b-a} \)
R1: \( (b-a) \)

Edited several times to correct english, missing sentences and code typos.

Encore une solution élégante de C.Ret
You play with LBL 1 through LBL 2 without RTN, subtle.
Tested on DM42. I forgot that there was a R^ on 15C.

[Gerson W. Barbosa]

HP-71B, HP-75B:

Code:


10 DATA 11,12,3,4,3,4,1,4,1,4
15 INPUT A,B
20 IF A<>B THEN Y=(A+B)/(A-B) ELSE P=2*PI*A @ GOTO 55
25 C=4/(2*Y-1)+1/3
30 FOR I=1 TO 5
35 READ N,D
40 C=N/(D*Y-C)
45 NEXT I
50 P=PI*(A-B)*(Y+C)
55 DISP P

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

where y = (a + b)/(a - b)

This approximation could be set up to return exact results for y = 1, but it would be needlessly longer.

I have yet to see how this approximation compares to Albert Chan's formula above. If I am not mistaken P(20, 1) and P(20, 0) in this approximation are better, but only very slightly. Also, every digit in his formula takes one step on the HP-15C, but it is probable that an HP-15C version of this program takes up even more.

The results on the HP-71B have been rounded to 10 figures in the table. An HP-15C version might return occasional differences in the last significant digit.


+-----+-----+-------+-------------+-------------+
|  a  |  b  |   h   |   result    |    exact    |
+-----+-----+-------+-------------+-------------+
| 20  | 20  | 0.000 | 125.6637061 | 125.6637061 |
| 20  | 19  | 0.026 | 122.5422527 | 122.5422527 |
| 20  | 18  | 0.053 | 119.4632087 | 119.4632087 |
| 20  | 17  | 0.081 | 116.4300496 | 116.4300496 |
| 20  | 16  | 0.111 | 113.4466716 | 113.4466716 | 
| 20  | 15  | 0.143 | 110.5174608 | 110.5174608 |
| 20  | 14  | 0.176 | 107.6473797 | 107.6473796 | 
| 20  | 13  | 0.212 | 104.8420720 | 104.8420720 |
| 20  | 12  | 0.250 | 102.1079955 | 102.1079955 | 
| 20  | 11  | 0.290 | 99.45258801 | 99.45258801 |
| 20  | 10  | 0.333 | 96.88448221 | 96.88448221 | 
| 20  |  9  | 0.379 | 94.41378489 | 94.41378488 |
| 20  |  8  | 0.429 | 92.05245041 | 92.05245038 | 
| 20  |  7  | 0.481 | 89.81479158 | 89.81479146 |
| 20  |  6  | 0.538 | 87.71820187 | 87.71820139 |
| 20  |  5  | 0.600 | 85.78421973 | 85.78421775 | 
| 20  |  4  | 0.667 | 84.04018649 | 84.04017816 |
| 20  |  3  | 0.739 | 82.52204269 | 82.52200588 | 
| 20  |  2  | 0.818 | 81.27965719 | 81.27948360 |
| 20  |  1  | 0.905 | 80.38937724 | 80.38851238 | 
| 20  |  0  | 1.000 | 80.03695848 | 79.99959924 | 
+-----+-----+-------|-------------+-------------+

———————

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)
…

[Albert Chan]

Here is Ellipse perimeter, AGM based version, without recursion.

10 DEF FNP(A,B) @ S=A*A+B*B @ T=1
20 REPEAT @ G=B @ K=(A-B)/2 @ B=SQRT(A*B) @ A=A-K @ T=T+T @ S=S-T*K*K @ UNTIL G=B
30 FNP = S/B*PI @ END DEF

> FNP(20,1)
80.3885123826


Comment: we can test convergence of S instead. This may reduce loop count by 1.

20 REPEAT @ K=(A-B)/2 @ B=SQRT(A*B) @ A=A-K @ T=T+T @ K=T*K*K @ S=S-K @ UNTIL S=S+K

[C.Ret]

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)  and R1:(a-b)
011-    RCL/1  2  x²  LSTx  GSB 0  3  1/x  +        ; Initiate  c_0_=4/(2*y- 1)+1/3 
019-    11  GSB.2  GSB 3  GSB 3  GSB 1  GSB 1       ; Iterate   c_i+1_=n/(d*y-c_i_)
026-    RCL+0  RCL*1  GTO 9                         ; Compute   (a-b)*(y+c)

Code:

029-     ►LBL.2    12  GTO 0                        ; ── ─── ─── d=12 ─┐    
033-     ►LBL 1     1  GTO 4                        ; ── n=1 ─┐        │
036-     ►LBL 3     3                               ; ── n=3 ─┤        │
038-     ►LBL 4     4                               ;         └─ d=4 ──┤
040-     ►LBL 0        RCL 0  *  R^  -  /  RTN      ;                  └─ c=n/(d*y-c) ───

Code:

047-  ►LBL 8     2  *                               ; Compute   P(a,a)=2πa
050-  ►LBL 9     PI  *                              ; Compute   P(a,b)~π(a-b)(y+1/(4y-1/(4y-3/(4y-3/(4y-11/(12y-(4/(2*Y-1)+1/3))))))) 
053- RTN