Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
12-10-2019, 11:50 PM (This post was last modified: 12-10-2019 11:55 PM by Gerson W. Barbosa.)
Post: #1
 Gerson W. Barbosa Senior Member Posts: 1,454 Joined: Dec 2013
Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
001:LBL  A ; Hypergeometric function –> F(a, b; c; z)
002:STOS  01
003:x⇆  Y
004:/
005:×
006:×
007:RCL  X
008:1
009:STO+  Y
010:x⇆  T
011:⇆  ZTYX
012:INC  04
013:INC  03
014:INC  02
015:INC  Y
016:RCL×  04
017:RCL/  02
018:RCL×  03
019:RCL/  Y
020:RCL×  01
021:RCL  Z
022:R↓
023:STO+  Z
024:⇆  TZXY
025:x≠?  Y
026:BACK  015
027:RTN

028:LBL  B ; Gauss-Kummer –> P = π(a + b)F(-1/2, -1/2; 1; h²)
029:©CONJ ; where h = (a - b)/(a + b)
031:x⇆  Y
033:STO  I
034:/
035:x²
036:#  1/2L
037:+/-
038:RCL  X
039:#  001
040:R↑
041:XEQ  A
042:RCL×  I
043:#  π
044:×  π
045:RTN

046:LBL  C ; Approximation formula using AGM
047:©ENTER; P ~ 2π{a + b - a*b/AGM(a,b) - 2[AGM(a,b) - √(a*b)]}
048:STO+  T
049:RCL×  Y
050:√
051:STO  I
052:x⇆  L
053:⇆  ZYXT
054:AGM
055:STO/  Y
056:RCL-  I
057:STO+  X
058:+
059:-
060:#  π
061:×
062:STO+  X
063:END

Examples:

8 ENTER 9 B -> 53.45328500297187553380768922447455

8 ENTER 9 C -> 53.45328500297(5636)

2 ENTER 3 B -> 15.86543958929058979133166302778306

2 ENTER 3 C -> 15.865439(6104)

Note: if a/b > 3, where a is the semi-major axis, then Cayley’s method (not implemented here) is more efficient. Please refer to this paper for more details:

http://web.tecnico.ulisboa.pt/~mcasquilh...llipse.pdf

——

A couple of hypergeometric function identities:

ln(1 + x) = xF(1, 1; 2; -x)

arcsin(x) = xF(1/2, 1/2; 3/2; x²)

Examples:

1 ENTER ENTER 2 ENTER 0.5 A 0.5 +/- * -> -0.69314718056
2 LN + -> 6e-34

0.5 ENTER ENTER 1.5 ENTER 0.25 A 0.5 * -> 0.5235987756
6 * π - -> 1e-32
12-11-2019, 01:47 AM
Post: #2
 Paul Dale Senior Member Posts: 1,750 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
Very nice. A function I didn't get to add to the 34S natively....

Pauli
12-13-2019, 04:08 AM (This post was last modified: 12-13-2019 04:11 AM by Gerson W. Barbosa.)
Post: #3
 Gerson W. Barbosa Senior Member Posts: 1,454 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
(12-11-2019 01:47 AM)Paul Dale Wrote:  Very nice. A function I didn't get to add to the 34S natively....

Thanks. It needn’t be that long, though.

001:LBL  A
002:STOS
003:x⇆  Y
004:/
005:×
006:×
007:#  001
008:x⇆  Y
009:RCL+  Y
010:y⇆  Z
011:x⇆  Y
012:x⇆  L
013:INC  04
014:INC  03
015:INC  02
016:INC  Z
017:RCL×  04
018:RCL/  02
019:RCL×  03
020:RCL/  Z
021:RCL×  01
022:RCL+  Y
023:x≠?  Y
024:BACK  013
025:RTN

Well, most likely it needn’t be this long either. Anyway, saving two steps inside the loop might make it even faster.

Gerson.
12-13-2019, 05:04 AM
Post: #4
 Paul Dale Senior Member Posts: 1,750 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
The STOS command takes an argument

Local registers might be useful to avoid corrupting the global ones.

Pauli
12-13-2019, 05:16 AM
Post: #5
 Gerson W. Barbosa Senior Member Posts: 1,454 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
(12-13-2019 05:04 AM)Paul Dale Wrote:  The STOS command takes an argument

I’ve been painstakingly copying and pasting from the iOS emulator one line at at time. All lines but the first one comes with garbage from previous lines I have to delete. Sometimes I undercorrect, sometimes I overcorrect...

002:STOS 01
01-11-2020, 01:20 AM (This post was last modified: 01-13-2020 02:59 AM by Gerson W. Barbosa.)
Post: #6
 Gerson W. Barbosa Senior Member Posts: 1,454 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
This new approximation formula involving AGM is more accurate than the previous one:

p ≈ 2π{agm(a,b)[77h⁴ - 768(h² - 1)] - [54h⁴ - 512(h² - 1)]√(ab)}/[23h⁴ - 256(h² - 1)]

where h = [(a - b)/(a + b)]

wp34s program:

001:LBL A
003:STO- Y
004:RCL+ T
005:/
006:x²
007:STO I
008:STO× I
009:DEC X
010:# 128
011:STO+ X
012:×
013:⇆ ZYXY
014:AGM
015:RCL L
016:RCL× T
017:√
018:# 054
019:RCL× I
020:RCL- T
021:RCL- T
022:×
023:# 077
024:RCL× I
025:RCL- T
026:RCL- T
027:RCL- T
028:RCL× Z
029:-
030:R↑
031:# 023
032:RCL× I
033:-
034:/
035:STO+ X
036:# π
037:×
038:RTN

Let’s try it with some examples:

9 ENTER 8 A -> 53.453285002971875533(62066)

9 ENTER 7 A -> 50.46202449950552(54160)

9 ENTER 6 A -> 47.596318767871(06954)

9 ENTER 4 A -> 42.36559853(28349)

9 ENTER 3 A -> 40.094679(4402)

9 ENTER 2 A -> 38.155(3995)

9 ENTER 1 A -> 36.68(6938)

Halley’s comet orbit:

a = 2667950000 km
b = 678281900 km

p ≈ 11464317964.011670 km
p = 11464318984.102995 km

difference: 1020.091 km

Let’s compare it with Ramanujan’s second approximation:

p ≈ π(a + b){1 + 3h²/[10 + √(4 - 3h²)]}

where h = [(a - b)/(a + b)]

039:LBL B
040:⇆ XYYZ
041:STO+ Z
042:-
043:RCL/ Y
044:x²
045:# 003
046:×
047:SDR 001
048:# 004
049:RCL- L
050:√
051:SDR 001
052:INC X
053:/
054:INC X
055:×
056:# π
057:×
058:END

9 ENTER 8 B -> 53.45328500297187(49239)

9 ENTER 7 B -> 50.46202449950(44277)

9 ENTER 6 B -> 47.596318767(75370)

9 ENTER 4 B -> 42.365598(45195)

9 ENTER 3 B -> 40.09467(8339)

9 ENTER 2 B -> 38.155(3915)

9 ENTER 1 B -> 36.687(5063)

Halley’s comet orbit:

a = 2667950000 km
b = 678281900 km

p ≈ 11464316399.11117 km
p = 11464318984.10300 km

difference: 2584.992 km

As the ellipse becomes more and more eccentric Ramanujan’s formula will eventually give smaller errors when compared to the AGM approximation, but overall the latter is better, at the cost of more complexity. For calculators with built-in AGM like the wp34s, that should be an option when Hypergeometric function is not available.

——-

PS: 01-13-2020, 02:59 AM

Here’s yet another approximation for the perimeter of the ellipse. It’s a bit less accurate than the previous one, but no square root is required, only basic arithmetic operations.

p ≈ π{(a + b)[h²(3h² - 136) + 320] + [a b/(a + b)](96h² - 256)}/[h²(3h² - 112) + 256]

where h = [(a - b)/(a + b)]

wp34s program:

001:LBL A
003:STO-Y
004:RCL+ T
005:/
006:x²
007:STO I
008:R↓
009:||
010:RCL L
011:RCL+ Z
012:# 003
013:RCL× I
014:# 136
015:-
016:RCL× I
017:# 160
018:+
019:RCL+ L
020:×
021:# 096
022:RCL× I
023:# 128
024:STO+ X
025:-
026:RCL× Z
027:+
028:# 003
029:RCL× I
030:# 112
031:-
032:RCL× I
033:# 128
034:+
035:RCL+ L
036:/
037:# π
038:×
039:END

9 ENTER 8 A -> 53.45328500297187(21835)

9 ENTER 7 A -> 50.462024499(4995106)

9 ENTER 6 A -> 47.596318767(23094)

9 ENTER 4 A -> 42.365598(09087)

9 ENTER 3 A -> 40.09467(3044)

9 ENTER 2 A -> 38.155(3243)

9 ENTER 1 A -> 36.68(6589)

Halley’s comet orbit:

a = 2667950000 km
b = 678281900 km

p ≈ 11464306668.94052 km
p = 11464318984.10299 km

difference: 12315.162 km
01-14-2020, 05:23 PM (This post was last modified: 01-15-2020 12:17 AM by Gerson W. Barbosa.)
Post: #7
 Gerson W. Barbosa Senior Member Posts: 1,454 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
Here is still another approximation to the perimeter of the ellipse:

p ≈ 2π[(aˢ + bˢ)/2]¹ᐟˢ

where

s = 3/2 + 1/(32/h² - ⅗h² - 4)

and where

h = [(a - b)/(a + b)]

That’s Thomas Muir’s approximation, I’ve only added an empirical correction term to his s = 3/2 parameter. The precision quite matches that of Ramanujan’s second approximation, the results being slightly above or below depending on the h parameter. In the case of the orbit of Halley’s comet, for instance, the absolute value of the difference is about three times less: 920.378 km. That’s about the flight distance from Warsaw to Stuttgart.

The wp34s program should be smaller (I’ll see to that later), but for convenience I will shift to Free42:

00 { 61-Byte Prgm }
01▸LBL "MB"
02 COMPLEX
03 ENTER
04 COMPLEX
05 ENTER
06 X<> ST Z
07 STO- ST Z
08 +
09 ÷
10 X↑2
11 32
12 X<>Y
13 ÷
14 0.6
15 RCL× ST L
16 -
17 4
18 -
19 1/X
20 1.5
21 +
22 X<>Y
23 COMPLEX
24 RCL ST Z
25 Y↑X
26 X<>Y
27 LASTX
28 Y↑X
29 +
30 2
31 ÷
32 X<>Y
33 1/X
34 Y↑X
35 STO+ ST X
36 PI
37 ×
38 END

Examples:

9 ENTER 8 XEQ MB ->

53.45328500297186507413258774441006

2 ENTER 3 XEQ MB ->

15.86543958923382742479398274914387

4 ENTER 1 XEQ MB ->

17.15684512619729903456675762188867

For exact results, use the ELP and 2F1 programs below, the latter copied and pasted from Jean-Marc Baillard in the Old HP-41C Software Library

00 { 38-Byte Prgm }
01▸LBL "ELP"
02 RCL- ST Y
03 X<>Y
04 RCL+ ST L
05 STO 04
06 ÷
07 X↑2
08 -0.5
09 RCL ST X
10 1
11 R↑
12 XEQ " 2F1"
13 RCL× 04
14 PI
15 ×
16 END

00 { 58-Byte Prgm }
01▸LBL " 2F1"
02 STO "Z"
03 R↓
04 STO 03
05 R↓
06 STO 02
07 R↓
08 STO 01
09 CLST
10 SIGN
11 ENTER
12 ENTER
13▸LBL 01
14 R↓
15 X<>Y
16 RCL 01
17 R↑
18 STO+ ST Y
19 R↓
20 ×
21 RCL 02
22 R↑
23 STO+ ST Y
24 R↓
25 ×
26 RCL 03
27 R↑
28 STO+ ST Y
29 ISG ST X
30 CLX
31 STO× ST Y
32 R↓
33 ÷
34 RCL× "Z"
35 STO ST Z
36 X<>Y
37 STO+ ST Y
38 X≠Y?
39 GTO 01
40 END

——-

PS: 01-15-2020, 00:17 AM

wp34s

001:LBL A
003:STO- Y
004:RCL+ T
005:/
006:x²
007:# 032
008:x⇆ Y
009:/
010:# 003
011:RCL× L
012:SDR 001
013:STO+ X
014:-
015:# 004
016:-
017:1/x
018:# 015
019:SDR 001
020:+
021:yᵡ
022:x⇆ Y
023:RCL L
024:yᵡ
025:STO+ Y
026:x⇆ L
027:# 002
028:STO/ Z
029:R↓
030:ᵡ√y
031:STO+ X
032:# π
033:×
034:END

If a = b, that is, when h = 0, the program will return “+∞ Error” (“Divide by 0” on the HP-42S). If this is an issue, just test for a = b in the beginning of the program and jump to current line # 31 if true. Or use this equivalent expression for s:

s = [h²(9h² + 50) - 480]/[h²(6h² + 40) - 320]
01-14-2020, 11:47 PM
Post: #8
 Albert Chan Senior Member Posts: 1,848 Joined: Jul 2018
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
(01-14-2020 05:23 PM)Gerson W. Barbosa Wrote:  Here is still another approximation to the perimeter of the ellipse:

p ≈ 2π[(aˢ + bˢ)/2]¹ᐟˢ

where

s = 3/2 + 1/(32/h² - ⅗h² - 4)

and where

h = [(a - b)/(a + b)]

How did you get "⅗" from the keyboard ?
BTW, if you change "⅗" to "217/360", above formula is always better than Ramanujan II

Let ellipse perimeter = Pi (a+b) corr(h)
To simplify, assume a=1, we have h = (1-b)/(1+b) → b = (1-h)/(1+h)

Below Mathematica code tried to compare error of corr(h):

In[1]:= corr[h_] := 4*EllipticE[1-b^2] / (Pi(1+b)) /. b->(1-h)/(1+h)
In[2]:= ramj2[h_]:= 1 + 3h^2/(10 + (4-3h^2)^(1/2))
In[3]:= Series[corr[h] - ramj2[h], {h,0,12}] // Simplify
﻿ ﻿ ﻿ ﻿ → (3/131072) h^10 + (79/2097152) h^12 + O(h^14)

In[4]:= muir[b_,s_] := 2*((1+b^s)/2)^(1/s) / (1+b)
In[5]:= muir2[h_]:= muir[(1-h)/(1+h), 3/2+1/(32/h^2 - 3/5*h^2 - 4)]
In[6]:= Series[corr[h] - muir2[h], {h,0,12}] // Simplify
﻿ ﻿ ﻿ ﻿ → (1/737280) h^8 + (41/13762560) h^10 + (-31319/825753600) h^12 + O(h^14)

In[7]:= muir3[h_]:= muir[(1-h)/(1+h), 3/2+1/(32/h^2 - 217/360*h^2 - 4)]
In[8]:= Series[corr[h] - muir3[h], {h,0,12}] // Simplify
﻿ ﻿ ﻿ ﻿ → (181/61931520) h^10 + (-16717/440401920) h^12 + O(h^14)

New constant removed h^8 term.
Perimeter error is smaller until h > 0.3946 (eccentricity > 0.9009).
01-15-2020, 12:59 AM (This post was last modified: 01-15-2020 12:59 AM by Gerson W. Barbosa.)
Post: #9
 Gerson W. Barbosa Senior Member Posts: 1,454 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
(01-14-2020 11:47 PM)Albert Chan Wrote:
(01-14-2020 05:23 PM)Gerson W. Barbosa Wrote:  Here is still another approximation to the perimeter of the ellipse:

p ≈ 2π[(aˢ + bˢ)/2]¹ᐟˢ

where

s = 3/2 + 1/(32/h² - ⅗h² - 4)

and where

h = [(a - b)/(a + b)]

How did you get "⅗" from the keyboard ?
BTW, if you change "⅗" to "217/360", above formula is always better than Ramanujan II

a) I simply copied it from Wikipedia ;

b) Great! Thanks for your analysis, results and experimenting with Mathematica.

I’ve used a much more modest tool, the hp33s Solver:

XROOT(X:(A^X+B^X)÷)=C

By making, for instance, A=3, B=2, C=2.52506313496 (the equivalent radius) and solving for X, I got X = 1.50125631949. The inverse of that result after the subraction of 3/2 is 795.975887056, which divided by 25 (1/h)² gives 31.8390348302. Likewise the same procedure on a couple of other examples always would return answers close to 32. Once this constant is settled, the second constant, 4, is evident: 25×32 - 795.975887056 = 4.02412294. And so on...
01-16-2020, 04:18 AM (This post was last modified: 01-16-2020 04:19 AM by Gerson W. Barbosa.)
Post: #10
 Gerson W. Barbosa Senior Member Posts: 1,454 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
(01-14-2020 11:47 PM)Albert Chan Wrote:  BTW, if you change "⅗" to "217/360", above formula is always better than Ramanujan II

Except at the vicinity of h = 1, when Ramanujan II becomes better, albeit only slightly.

Here are the fixed expression for s and the new HP-42S/Free42 version:

s = 3/2 + h²/[32 - h²(217h²/360 + 4)]

00 { 66-Byte Prgm }
01▸LBL "MB"
02 COMPLEX
03 ENTER
04 COMPLEX
05 ENTER
06 X<> ST Z
07 STO- ST Z
08 +
09 ÷
10 X↑2
11 0.361
12 →HR
13 RCL× ST Y
14 4
15 +
16 RCL× ST Y
17 +/-
18 32
19 +
20 ÷
21 1.5
22 +
23 X<>Y
24 COMPLEX
25 RCL ST Z
26 Y↑X
27 X<>Y
28 LASTX
29 Y↑X
30 +
31 2
32 ÷
33 X<>Y
34 1/X
35 Y↑X
36 STO+ ST X
37 PI
38 ×
39 END

The error when computing the length of Halley’s comet is now 1146.924 km, about the flight distance from Milan to Warsaw („Z ziemi włoskiej do Polski”), but the results are overall better.
01-18-2020, 01:14 AM (This post was last modified: 01-18-2020 01:19 AM by Gerson W. Barbosa.)
Post: #11
 Gerson W. Barbosa Senior Member Posts: 1,454 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
Here is another approximation that relies only on basic arithmetic operations and takes up less than 40 steps on HP-42S/Free42, not counting the final END:

p ≈ π(a + b){h²[h²(-93h² + 224) + 2304] - 4096}/{h²[h²(7h² - 544) + 3328] - 4096}

00 { 80-Byte Prgm }
01▸LBL "EL"
02 RCL+ ST Y
03 STO 01
04 X<>Y
05 LASTX
06 RCL- ST Y
07 X<>Y
08 RCL+ ST L
09 ÷
10 X↑2
11 ENTER
12 ENTER
13 ENTER
14 -93
15 ×
16 224
17 +
18 ×
19 2304
20 +
21 ×
22 4096
23 -
24 STO× 01
25 R↓
26 7
27 ×
28 544
29 -
30 ×
31 3328
32 +
33 ×
34 4096
35 -
36 X<> 01
37 RCL÷ 01
38 PI
39 ×
40 END

The error in the usual example (orbit of Halley’s comet) is 477.529 km, the flight distance from Rome to Milan.
Use the formula or the program to check the error in the orbit of former planet Pluto (a = 5906376272 km, b = 5720637952.8 km).
01-18-2020, 10:25 AM
Post: #12
 grsbanks Senior Member Posts: 1,219 Joined: Jan 2017
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
Could you enclose your listings in "code" tags? It makes the whole post shorter and tidier. See below:

(01-18-2020 01:14 AM)Gerson W. Barbosa Wrote:
Code:
00 { 80-Byte Prgm } 01▸LBL "EL" 02 RCL+ ST Y 03 STO 01 04 X<>Y 05 LASTX 06 RCL- ST Y 07 X<>Y 08 RCL+ ST L 09 ÷ 10 X↑2 11 ENTER 12 ENTER 13 ENTER 14 -93 15 × 16 224 17 + 18 × 19 2304 20 + 21 × 22 4096 23 - 24 STO× 01 25 R↓ 26 7 27 × 28 544 29 - 30 × 31 3328 32 + 33 × 34 4096 35 - 36 X<> 01 37 RCL÷ 01 38 PI 39 × 40 END

There are only 10 types of people in this world. Those who understand binary and those who don't.
01-18-2020, 03:41 PM (This post was last modified: 01-19-2020 02:09 AM by Gerson W. Barbosa.)
Post: #13
 Gerson W. Barbosa Senior Member Posts: 1,454 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
(01-18-2020 10:25 AM)grsbanks Wrote:  Could you enclose your listings in "code" tags? It makes the whole post shorter and tidier.

As you please, but you’ll have to scroll down to see the new HP-42S/Free42 version is now 36 steps long, LBL and END included :-)

Code:
 00 { 73-Byte Prgm } 01▸LBL "EL" 02 RCL- ST Y 03 X<>Y 04 RCL+ ST L 05 STO 01 06 ÷ 07 X↑2 08 ENTER 09 ENTER 10 ENTER 11 7 12 × 13 544 14 - 15 × 16 3328 17 + 18 × 19 4096 20 - 21 STO÷ 01 22 R↓ 23 -93 24 × 25 224 26 + 27 × 28 2304 29 + 30 × 31 4096 32 - 33 RCL× 01 34 PI 35 × 36 END

The corresponding wp34s version is 41 steps longs – too many constants, remove a few and it will be perfect! :-)
Perhaps someone has a trick to save a step or two...

wp34s version:

Code:
 001:LBL A 002:⇆ XYYY 003:STO+ Y 004:©- 005:RCL/ I 006:x² 007:FILL 008:# 007 009:× 010:# 032 011:STO J 012:# 017 013:× 014:-   015:×   016:# 104 017:RCL× J 018:+ 019:× 020:# 128 021:RCL× J 022:- 023:STO/ I 024:R↓ 025:# 007 026:RCL× J 027:# 093 028:RCL× Z 029:- 030:× 031:# 072 032:RCL× J 033:+ 034:× 035:# 128 036:RCL× J 037:- 038:RCL× I 039:# π  040:× 041:END

p ≈ π(a + b){h²[h²(-93h² + 224) + 2304] - 4096}/{h²[h²(7h² - 544) + 3328] - 4096}

where

h = [(a - b)/(a + b)]

Thanks for your suggestion. This indeed looks cleaner now.

Edited after saving one step in the wp34s program.

PS:

1) The wp34s program can fit in 40 steps, including LBL and END;

2) The following is even more accurate, giving an error of only 116.385 km for our usual example – that’s about the distance from Milan to Parma.

p ≈ π(a + b){h²[h²(-381h² + 548) + 8448] - 13312}/{h²[h²(34h² - 2188) + 11776] - 13312}

Same number of steps on Free42/HP-42S, taking up a few more bytes.

Code:
 00 { 79-Byte Prgm } 01▸LBL "EL" 02 RCL- ST Y 03 X<>Y 04 RCL+ ST L 05 STO 01 06 ÷ 07 X↑2 08 ENTER 09 ENTER 10 ENTER 11 34 12 × 13 2188 14 - 15 × 16 11776 17 + 18 × 19 13312 20 - 21 STO÷ 01 22 R↓ 23 -381 24 × 25 548 26 + 27 × 28 8448 29 + 30 × 31 13312 32 - 33 RCL× 01 34 PI 35 × 36 END
01-21-2020, 09:01 PM
Post: #14
 Gerson W. Barbosa Senior Member Posts: 1,454 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
(01-18-2020 03:41 PM)Gerson W. Barbosa Wrote:  p ≈ π(a + b){h²[h²(-93h² + 224) + 2304] - 4096}/{h²[h²(7h² - 544) + 3328] - 4096}

where

h = [(a - b)/(a + b)]

...

1) The wp34s program can fit in 40 steps, including LBL and END;

Actually it does fit in 38 steps:

Code:
 001:LBL A 002:⇆ XYYY 003:STO+ Y 004:©- 005:RCL/ I 006:x² 007:x³ 008:STO J 009:x⇆ L 010:FILL 011:# 007 012:×  013:# 072 014:+  015:XEQ 00 016:# 093 017:RCL× J 018:- 019:RCL× I 020:# 104 021:# 017 022:RCL× T 023:- 024:XEQ 00 025:# 007 026:RCL× J 027:+ 028:/ 029:# π   030:×  031:RTN 032:LBL 00 033:RCL× Z 034:# 128 035:-  036:# 032 037:×  038:END

Since this is the only wp34s program which is more or less optimized, I will use it as the basis for the AGM method present by Albert Chan here:

Code:
 001:LBL A 002:RCL× Y 003:z⇆ L  004:# π 005:x⇆ Y  006:× 007:STO 01 008:x⇆ L 009:√ 010:⇆ ZYYX 011:AGM 012:STO/ 01 013:x⇆ L 014:+ 015:2 016:/ 017:⇆ XYYY 018:STO+ Y 019:©- 020:RCL/ I 021:x² 022:x³ 023:STO J 024:x⇆ L 025:FILL 026:# 007 027:×  028:# 072 029:+  030:XEQ 00 031:# 093 032:RCL× J 033:- 034:RCL× I 035:# 104 036:# 017 037:RCL× T 038:- 039:XEQ 00 040:# 007 041:RCL× J 042:+ 043:/ 044:# π   045:×  046:RCL- 01 047:STO+ X 048:RTN 049:LBL 00 050:RCL× Z 051:# 128 052:-  053:# 032 054:×  055:END

Example:

Halley’s comet orbit

a = 2667950000 km

b = 678281900 km

2667950000 ENTER 678281900 A ->

p ≈ 11464318984.10299492347510528221642 km

p = 11464318984.10299540114254189105734 km

difference: 477.667 µm
01-25-2020, 02:53 AM
Post: #15
 Gerson W. Barbosa Senior Member Posts: 1,454 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
(01-21-2020 09:01 PM)Gerson W. Barbosa Wrote:   difference: 477.667 µm

Not even NASA needs to compute the perimeter of an ellipse the shape and dimension of the orbit of Halley’s comet to a precision of less than half a millimeter.
So, in order to minimize program size I’ve worked out another approximation which will give reasonable results when combined with the AGM technique, despite its sloppiness.

p ≈ π(a + b)[(h² + 8)/8]²

where

h = [(a - b)/(a + b)]

Code:
 001:LBL A 002:©ENTER 003:× 004:STO I 005:√ 006:x⇆ T 007:RCL+ Y 008:R↓ 009:AGM 010:STO/ I 011:x⇆ Z  012:2 013:/ 014:RCL+ Y 015:RCL L 016:RCL- Z 017:RCL/ Y 018:x² 019:# 008 020:+  021:RCL/ L 022:x² 023:× 024:# π  025:STO× I 026:× 027:RCL- I 028:STO+ X 029:END

The difference in our usual example is now 122.324 km and the error in the perimeter of the orbit of Pluto is 19.256 fm (femtometer), still beyond any practical requirement. Anyway, the program fits in only 29 steps.

The longer version is now 53 steps long:

Code:
 001:LBL A 002:©ENTER 003:× 004:STO 01 005:√ 006:x⇆ T 007:RCL+ Y 008:R↓ 009:AGM 010:STO/ 01 011:x⇆ Z 012:2 013:/ 014:⇆ XYYY 015:STO+ Y 016:©- 017:RCL/ I 018:x² 019:x³ 020:STO J 021:x⇆ L 022:FILL 023:# 007 024:×  025:# 072 026:+  027:XEQ 00 028:# 093 029:RCL× J 030:- 031:RCL× I 032:# 104 033:# 017 034:RCL× T 035:- 036:XEQ 00 037:# 007 038:RCL× J 039:+ 040:/ 041:# π   042:STO× 01 043:× 044:RCL- 01 045:STO+ X 046:RTN 047:LBL 00 048:RCL× Z 049:# 128 050:-  051:# 032 052:×  053:END
01-29-2020, 11:44 AM
Post: #16
 Ángel Martin Senior Member Posts: 1,324 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
(12-11-2019 01:47 AM)Paul Dale Wrote:  Very nice. A function I didn't get to add to the 34S natively....

I "discovered" it a few years back on JM Baillard's pages and have used a MCODE version of it it profusely in the SandMath module.

The applicability of this function is amazing, take a look at JM's nice documentation here: http://hp41programs.yolasite.com/fracintegrdiff.php
01-29-2020, 12:29 PM
Post: #17
 Paul Dale Senior Member Posts: 1,750 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
Yes, it covers a wide range of other functions.

The Jacobi elliptical functions were another suite I'd have liked to include.

Pauli
02-08-2020, 11:17 PM (This post was last modified: 02-09-2020 01:37 AM by Gerson W. Barbosa.)
Post: #18
 Gerson W. Barbosa Senior Member Posts: 1,454 Joined: Dec 2013
RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)
Here is an explanation of the (unorthodox) method I've been using to find these approximations. It is based on sheer observation rather than on complicated mathematical analysis. For example, The ratio between four times the equivalent radius and its difference from the arithmetic mean of the semi-axes is about sixteen times the inverse of the square of the h parameter (difference of the semi-axes over their sum) plus 3:

$$\frac{4r}{r-\frac{a+b}{2}}\approx \frac{16}{h^{2}}+3$$

By solving it for r we obtain

$$r \approx \frac{a+b}{2}\left ( \frac{3h^{2}+16}{16-h^{2}} \right )$$

or

$$p\approx \pi \left ( a+ b \right )\left ( \frac{3h^{2}+16}{16-h^{2}} \right )$$

(1 655 947.321 km) (Numbers attached to a length unity between parentheses here and elsewhere following a formula or program refer to the error obtained when using it to compute the perimeter of the orbit of Halley's comet).

That's Selmer's approximation from 1975 (Gauss, Landen, Ramanujan, the Arithmetic-Geometric Mean, Ellipses, π, and the Ladies Diary, page 600). This linked paper contains a derivation of Ramanujan's second approximation (page 602).

This can be carried one step further:

$$\frac{4r }{r - \frac{a+b}{2}} \approx c - \frac{1}{\frac{c-9}{3}}$$

where

$$c = \frac{16}{h^{2}}+3$$

$$r \approx \frac{a+b}{2} \left ( \frac{256-48h^{2}-21h^{4}}{256-112h^{2}+3h^{4} } \right )$$

or

$$p \approx \pi \left ( a+b \right )\left ( \frac{256-48h^{2}-21h^{4}}{256-112h^{2}+3h^{4} } \right )$$

(12 315.162 km)

That's Jacobsen's and Waadeland's approximation from 1985 (page 601 in the previously linked paper).

That can be carried even further:

$$\frac{4r }{r - \frac{a+b}{2}} \approx c - \frac{1}{\frac{c-9}{3}-\frac{1}{\frac{9d}{11}}}$$

where

$$c = \frac{16}{h^{2}}+3$$

and

$$d = \frac{\left ( \frac{16}{h^{2}}+3 \right )-9}{3 }$$

which implies

$$r \approx\left ( \frac{a+b}{2} \right )\left ( \frac{4096-2304h^{2}-224h^{4}+93h^{6}}{4096-3328h^{2}+544h^{4}-7h^{6}} \right )$$

or

$$p \approx \pi \left ( a+b \right )\left ( \frac{4096-2304h^{2}-224h^{4}+93h^{6}}{4096-3328h^{2}+544h^{4}-7h^{6}} \right )$$

(477.529 km)

-----------------

The next three formulae are not particularly useful, but they might be interesting nonetheless:

-----------------

$$p\approx \pi \left ( a+b \right )\left ( 1+\frac{1}{8h^{-2}-\frac{1}{8h^{-2}-\frac{17}{8}}} \right )^{2}$$

(23 933.122 km)

-----------------

$$p\approx 2\pi \left [ a+b- AGM\left (a-c,b+c \right ) \right ]$$

where

$$c=\left ( a-b \right )\left ( \frac{h^{2}}{16}+\frac{h^{4}}{256}+\frac{h^{6}}{4096} \right )$$

$$h=\frac{a-b}{a+b }$$

and

$$a\geqslant b$$

Code:
 0001 **LBL A 0002 x<? Y 0003 x[<->] Y 0004 [cmplx]STO A 0005 [<->] XYYZ 0006 STO+ Z 0007 - 0008 [cmplx]ENTER 0009 RCL/ Y 0010 x[^2] 0011 # 016 0012 / 0013 ENTER[^] 0014 INC X 0015 RCL[times] Y 0016 INC X 0017 [times] 0018 RCL[times] Z 0019 RCL+ A 0020 RCL B 0021 RCL- L 0022 AGM 0023 - 0024 STO+ X 0025 # [pi] 0026 [times] 0027 END

(8040.370 km)

-----------------

$$p\approx 2\pi\left ( \frac{\frac{\left ( a-b \right )^{2}}{4h\cdot AGM\left ( a,b \right )}}{1+\frac{h^{4}}{8-4h^{2}-\frac{25h^{4}}{16}-\frac{h^{6}}{\frac{32}{9}-\frac{9h^{2}}{4}}}} \right )$$

(4754.164 km)

-----------------

It appears the best option in terms of size, speed and accuracy is the AGM method presented by Albert Chan associated with Ramanujan's second approximation, which requires no more than 35 steps on the wp34s:

Code:
 0001 **LBL A 0002 [cmplx]ENTER 0003 [times] 0004 STO I 0005 [sqrt] 0006 x[<->] T 0007 RCL+ Y 0008 R[v] 0009 AGM 0010 STO/ I 0011 x[<->] Z 0012 # 002 0013 / 0014 [<->] XYYZ 0015 STO+ Z 0016 - 0017 RCL/ Y 0018 x[^2] 0019 # 003 0020 [times] 0021 SDR 001 0022 # 004 0023 RCL- L 0024 [sqrt] 0025 SDR 001 0026 INC X 0027 / 0028 INC X 0029 [times] 0030 # [pi] 0031 STO[times] I 0032 [times] 0033 RCL- I 0034 STO+ X 0035 END

(10.132 cm)