Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s)

[Gerson W. Barbosa]

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)
030:©RCL  L
031:x⇆  Y 
032:©+
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

[Paul Dale]

Very nice. A function I didn't get to add to the 34S natively....

Pauli

[Gerson W. Barbosa]

(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.

[Paul Dale]

The STOS command takes an argument

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

Pauli

[Gerson W. Barbosa]

(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...
It should read

002:STOS 01

[Gerson W. Barbosa]

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
002:©ENTER
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
002:©ENTER
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

[Gerson W. Barbosa]

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
002:©ENTER
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]

[Albert Chan]

(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).

[Gerson W. Barbosa]

(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...

[Gerson W. Barbosa]

(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.

[Gerson W. Barbosa]

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

[grsbanks]

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

[Gerson W. Barbosa]

(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

[Gerson W. Barbosa]

(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

[Gerson W. Barbosa]

(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

[Ángel Martin]

(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

[Paul Dale]

Yes, it covers a wide range of other functions.

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


Pauli

[Gerson W. Barbosa]

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

which leads to

\(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)

Edited to fix a hyperlink