+- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: Not HP Calculators (/forum-7.html) +--- Forum: Not quite HP Calculators - but related (/forum-8.html) +--- Thread: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) (/thread-14148.html) Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Gerson W. Barbosa - 12-10-2019 11:50 PM 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/~mcasquilho/compute/com/,ellips/PerimeterOfEllipse.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 RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Paul Dale - 12-11-2019 01:47 AM Very nice. A function I didn't get to add to the 34S natively.... Pauli RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Gerson W. Barbosa - 12-13-2019 04:08 AM (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. RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Paul Dale - 12-13-2019 05:04 AM The STOS command takes an argument Local registers might be useful to avoid corrupting the global ones. Pauli RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Gerson W. Barbosa - 12-13-2019 05:16 AM (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 RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Gerson W. Barbosa - 01-11-2020 01:20 AM 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 RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Gerson W. Barbosa - 01-14-2020 05:23 PM 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] RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Albert Chan - 01-14-2020 11:47 PM (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). RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Gerson W. Barbosa - 01-15-2020 12:59 AM (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... RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Gerson W. Barbosa - 01-16-2020 04:18 AM (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. RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Gerson W. Barbosa - 01-18-2020 01:14 AM 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). RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - grsbanks - 01-18-2020 10:25 AM 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 RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Gerson W. Barbosa - 01-18-2020 03:41 PM (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 RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Gerson W. Barbosa - 01-21-2020 09:01 PM (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 RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Gerson W. Barbosa - 01-25-2020 02:53 AM (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 RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Ángel Martin - 01-29-2020 11:44 AM (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 RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Paul Dale - 01-29-2020 12:29 PM Yes, it covers a wide range of other functions. The Jacobi elliptical functions were another suite I'd have liked to include. Pauli RE: Hypergeometric function – Perimeter of an Ellipse and other applications (wp34s) - Gerson W. Barbosa - 02-08-2020 11:17 PM 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 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