+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (28/48/50) Complete Elliptic Integrals (/thread-21364.html)
(28/48/50) Complete Elliptic Integrals - John Keith - 02-26-2024 08:14 PM
Note: All of the programs have been changed. E(k) has been removed and two new programs added.
These programs compute the complete elliptic integrals of the first kind K(k), second kind E(k) and third kind Π(n, k). Also the exact perimeter of the ellipse with semi-major axis a and semi-minor axis b, and the elliptic nome. They are reasonably fast (less than 300ms on my 50g) and accurate with 11 or 12 correct digits for almost all inputs, real or complex.
The first program Kk returns the complete elliptic integral of the first kind. The second program EKk returns E(k) on level 2 and K(k) on level 1, similar to the Matlab function ellipke. The reason for this program is that the computation of Ek does most of the work of computing Kk, and this program is significantly faster than the other two separately. The third program Π(n, k) computes the complete elliptic integral of the third kind given n on level 2 and k on level 1. It is a rather large program due to checks for several special cases. Note: for the HP-28, # 203h DOERR should be replaced by a string, e.g. "Bad Argument Value", since the HP-28 lacks the DOERR command.
The next program ELPER returns the perimeter of the ellipse given a on level 2 and b on level 1. If a and b are reversed, the result is a complex number with the perimeter as the real part and 0 as the imaginary part. The last program ENOME computes the elliptic nome q(x), a function related to the elliptic integral of the first kind.
These programs are based on formulae in this Wikipedia page, and the DLMF. Hugh Steers' C program in this thread from the old forum computes the perimeter in a slightly different way but seems to be essentially the same algorithm. Kk and EKk have been improved based on Albert Chan's suggestion in post #7 below.
Important note: the equivalent Matlab and Mathematica functions compute K(m) and E(m), where m = k^2, so the inputs need to be adjusted accordingly when comparing results.
The programs are listed in the form of a directory due to several dependencies.
Code:
DIR
Kk
\<< DUP ABS 1
IF \=/
THEN 1 1 ROT SQ - \v/
DO DUP ROT ROT DUP2 + 2 / ROT ROT * \v/
UNTIL ROT OVER ==
END DROP 2 * \pi \->NUM SWAP /
ELSE DROP MAXR \->NUM
END
\>>
EKk
\<< DUP ABS 1
IF \=/
THEN 1 SWAP SQ DUP2 - \v/
ROT 3 PICK 2 / 4 ROLLD .5 ROT ROT
WHILE ROT 2 * ROT ROT
DUP ROT ROT DUP2 * \v/
ROT ROT + 2 / ROT OVER \=/
REPEAT 4 ROLL OVER 4 * / SQ
4 PICK OVER * 6 ROLL +
5 ROLLD 4 ROLLD
END 5 ROLLD DROP2 DROP
\pi \->NUM SWAP 1 SWAP - OVER *
3 PICK 2 * / ROT 2 * ROT SWAP /
ELSE DROP 1 MAXR \->NUM
END
\>>
\PInk
\<<
IF OVER ABS
THEN
IF OVER 1 == OVER ABS 1 == OR
THEN DROP2 MAXR \->NUM
ELSE
IF OVER DUP RE 1 > SWAP IM NOT AND
THEN DROP2 # 203h DOERR
ELSE SWAP 1 1 \-> n s q
\<< 1 1 ROT SQ - \v/ OVER n - ROT ROT
DO 3 DUPN * DUP2 + ROT ROT - OVER /
q * 2 / s 6 ROLLD 6 PICK
OVER + 's' STO 'q' STO
4 ROLL \v/ 2 * / SQ ROT ROT
DUP ROT ROT DUP2 + 2 / ROT ROT * \v/
UNTIL ROT OVER == 5 ROLL s == AND
END ROT ROT DROP2
s n 1 OVER - / * 2 +
\pi \->NUM * SWAP 4 * /
\>>
END
END
ELSE NIP Kk
END
\>>
ELPER
\<< OVER / SQ 1 SWAP - \v/ EKk DROP * 4 *
\>>
ENOME
\<< 1 OVER SQ - \v/ Kk \pi NEG \->NUM * SWAP Kk / EXP
\>>
ENDAnd finally a listing optimized for the HP 49 and 50 with decimal points and extra stack commands.
Code:
DIR
Kk
\<< DUP ABS 1.
IF \=/
THEN 1. 1. ROT SQ - \v/
DO DUP UNROT DUP2 + 2. /
UNROT * \v/
UNTIL ROT OVER ==
END DROP 2. * \pi \->NUM SWAP /
ELSE DROP MAXR \->NUM
END
\>>
EKk
\<< DUP ABS 1.
IF \=/
THEN 1. SWAP SQ DUP2 - \v/
ROT PICK3 2. / 4. ROLLD .5 UNROT
WHILE ROT 2. * UNROT
DUP UNROT DUP2 * \v/
UNROT + 2. /
ROT OVER \=/
REPEAT 4. ROLL OVER 4. * / SQ
4. PICK OVER * 6. ROLL +
5. ROLLD 4. ROLLD
END 5. ROLLD DROP2 DROP
\pi \->NUM SWAP 1. SWAP - OVER *
PICK3 2. * / ROT 2. * ROT SWAP /
ELSE DROP 1. MAXR \->NUM
END
\>>
\PInk
\<< IF OVER ABS
THEN
IF OVER 1. == OVER ABS 1. == OR
THEN DROP2 MAXR \->NUM
ELSE
IF OVER DUP RE 1. > SWAP IM NOT AND
THEN DROP2 #203h DOERR
ELSE SWAP 1. 1. \-> n s q
\<< 1. 1. ROT SQ - \v/ OVER n - UNROT
DO PICK3 PICK3 PICK3 *
DUP2 + UNROT - OVER /
q * 2. / s 6. ROLLD 6. PICK
OVER + 's' STO 'q' STO
4. ROLL \v/ 2. * / SQ UNROT
DUP UNROT DUP2 + 2. / UNROT * \v/
UNTIL ROT OVER == 5. ROLL s == AND
END UNROT DROP2
s n 1. OVER - / * 2. +
\pi \->NUM * SWAP 4. * /
\>>
END
END
ELSE NIP Kk
END
\>>
ELPER
\<< OVER / SQ 1. SWAP - \v/ EKk DROP * 4. *
\>>
ENOME
\<< 1. OVER SQ - \v/ Kk \pi NEG \->NUM * SWAP Kk / EXP
\>>
ENDRE: (28/48/50) Complete Elliptic Integrals - Thomas Klemm - 02-27-2024 12:08 AM
Just out of curiosity: why would you use .00000000002 instead of 2E-11?
Is this faster?
From Legendre's relation:
For example:
\(
K({\color {blueviolet}{\tfrac {3}{5}}})E({\color {blue}{\tfrac {4}{5}}})+E({\color {blueviolet}{\tfrac {3}{5}}})K({\color {blue}{\tfrac {4}{5}}})-K({\color {blueviolet}{\tfrac {3}{5}}})K({\color {blue}{\tfrac {4}{5}}})={\tfrac {1}{2}}\pi
\)
.6 Kk
.8 Ek
*
.6 Ek
.8 Kk
*
+
.6 Kk
.8 Kk
*
-
1.57079632679
Compare this to \(\frac{\pi}{2}\):
1.5707963268
LGTM
RE: (28/48/50) Complete Elliptic Integrals - John Keith - 02-27-2024 08:53 PM
(02-27-2024 12:08 AM)Thomas Klemm Wrote: Just out of curiosity: why would you use .00000000002 instead of 2E-11?
Is this faster?
It's exactly the same; if you type 2E-11 into the calculator, it displays .00000000002. I suppose it would be better to have 2E-11 in the listings above to make life easier for human readers.
RE: (28/48/50) Complete Elliptic Integrals - Gjermund Skailand - 02-28-2024 11:21 AM
There is also my old post in the old forum which you may find interesting
: http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv021.cgi?read=245627
br Gjermund
RE: (28/48/50) Complete Elliptic Integrals - floppy - 02-28-2024 03:57 PM
Using AGM and MAGM is perhaps quicker?
https://stackoverflow.com/questions/42310956/how-to-calculate-the-perimeter-of-an-ellipse/69606140#69606140
RE: (28/48/50) Complete Elliptic Integrals - John Keith - 02-28-2024 09:54 PM
(02-28-2024 03:57 PM)floppy Wrote: Using AGM and MAGM is perhaps quicker?
https://stackoverflow.com/questions/42310956/how-to-calculate-the-perimeter-of-an-ellipse/69606140#69606140
For the perimeter of the ellipse perhaps, but these programs compute the elliptic integrals which may be used for other purposes.
RE: (28/48/50) Complete Elliptic Integrals - Albert Chan - 02-29-2024 04:31 AM
Hi, John Keith
It may be better not to test |AM-GM| below some epsilon. (2E-11)
Instead, test for equality of convergence of AM (or GM), see Werner's AGM code
This make AGM(a,b) code work for any sized (a,b), any machine precisions.
Modify the AGM, we can get both K and E at the same time.
Note: below Free42 "K" code, argument is parameter m = k^2
(06-20-2020 04:23 PM)Albert Chan Wrote: x, y = agm2(a, b)
→ x = converged GM of agm(a, b)
→ y = -Σ(2k (½gapk)² , k = 1 .. n), n = number of iterations to converge GM
With this new setup, ellipse_perimeter(a,b) = 4 a E(1-(b/a)²) = pi (y + b² + a²)/x
Example: for ellipse perimeter, a=50, b=10
50 Enter 10 XEQ "AGM"
[X<>Y] 10 [X↑2] + 50 [X↑2] + [X<>Y] ÷ ; ellipse "diameter" ≈ 66.8770488614
PI × ; ellipse perimeter ≈ 210.100445397
For E(m), K(m):
x, y = agm2(1, sqrt(1-m))
→ K = pi / (2*x)
→ E = K * (y + (1-m) + 1²) / 2 = K + K*(y-m)/2
Instead of returning (E, K), it is better to return (K, E-K), keeping the accurate difference.
(for 0 < m < 1, y is negative, E-K = K*(y-m)/2 avoided subtraction cancellation errors)
Again, for ellipse_perimeter, a=50, b=10, e² = 1-(b/a)² = 0.96
0.96 XEQ "K" + 200 × ; ellipse perimeter ≈ 210.100445397
Code:
00 { 92-Byte Prgm }
01▸LBL "K" ; (m) -> K, E-K, m, m
02 LSTO "m"
03 1
04 ENTER
05 RCL- "m" ; 1-m 1
06 SQRT
07 XEQ "AGM" ; x y
08 PI
09 X<>Y
10 STO+ ST X ; x+x pi y
11 ÷
12 RCL "m" ; m K y
13 STO- ST Z
14 X<> ST Z ; y-m K m
15 2
16 ÷
17 X<>Y
18 STO× ST Y ; K E-K m m
19 RTN
20▸LBL "AGM" ; (a,b) -> (agm, acc)
21 LSTO "A"
22 CLX
23 LSTO "S" ; s=0
24 SIGN
25 LSTO "T" ; t=1
26 X<>Y
27▸LBL 01 ; b .
28 ENTER
29 RCL× "A" ; ab b
30 LASTX
31 RCL- "A" ; b-a ab b
32 2
33 STO× "T"
34 ÷ ; k ab b b
35 STO+ "A"
36 X↑2
37 RCL× "T"
38 STO- "S"
39 R↓
40 SQRT ; GM b b tkk
41 X≠Y?
42 GTO 01
43 RCL "S"
44 X<>Y ; GM S b b
45 END
RE: (28/48/50) Complete Elliptic Integrals - John Keith - 02-29-2024 09:23 PM
Thanks, Albert. I will try your ideas and Werner's, although it will complicate the program- RPL has neither GOTO nor BREAK.
RE: (28/48/50) Complete Elliptic Integrals - Thomas Klemm - 02-29-2024 11:43 PM
(02-29-2024 09:23 PM)John Keith Wrote: RPL has neither GOTO nor BREAK.
Cf. RPL equivalents to BREAK, CYCLE, EXIT, etc.?
RE: (28/48/50) Complete Elliptic Integrals - John Keith - 03-04-2024 05:37 PM
Thanks, Thomas, I am aware of those methods, I was really just grumbling about having to restructure the program. It turned out to be fairly simple, just changing the DO loop to a WHILE loop.
RE: (28/48/50) Complete Elliptic Integrals - John Keith - 03-04-2024 06:11 PM
(02-29-2024 04:31 AM)Albert Chan Wrote: Modify the AGM, we can get both K and E at the same time.
Note: below Free42 "K" code, argument is parameter m = k^2
(06-20-2020 04:23 PM)Albert Chan Wrote: For E(m), K(m):
x, y = agm2(1, sqrt(1-m))
→ K = pi / (2*x)
→ E = K * (y + (1-m) + 1²) / 2 = K + K*(y-m)/2
Instead of returning (E, K), it is better to return (K, E-K), keeping the accurate difference.
(for 0 < m < 1, y is negative, E-K = K*(y-m)/2 avoided subtraction cancellation errors)
I tried your agm2 code and it gives the correct value for K(k), but does not seem to give the correct value for E(k) using your formula above.
Here is my RPL translation of your agm2:
Code:
\<< IF DUP ABS @ b != 0?
THEN OVER SQ UNROT -1. UNROT @ s, t
WHILE ROT 2. * UNROT DUP2 + 2. / @ t = 2*t, a1
PICK3 OVER \=/ @ a != a1?
REPEAT PICK3 ROT * \v/ @ New b
OVER 4. ROLL - SQ @ k^2
4. PICK * 5. ROLL + 4. ROLLD @ Update s
END DROP 4. ROLLD DROP2 @ Return b, s
ELSE SWAP SQ SWAP @ Else return a^2, 0
END
\>>For m = .25, I get agm2(1, sqrt(.75)) ~ 0.866025403784 which returns
b = 0.93180839162
s = 0.991019606076
For K(m) I get 1.68575035482 which is correct in all 12 digits and agrees with my program for K(.5) given that m = k^2. However, using your formula E = K + K*(y-m)/2 I get 2.31033738676, whereas the correct value is 1.46746220934.
It is entirely possible that I made a transcription error and the above program is returning the wrong value for s, but I can't see where the problem is.
Edit: Problem solved, see post below. Program edited to fix logical error.
RE: (28/48/50) Complete Elliptic Integrals - Albert Chan - 03-04-2024 06:39 PM
(03-04-2024 06:11 PM)John Keith Wrote: For m = .25, I get agm2(1, sqrt(.75) ~ 0.866025403784) which returns
b = 0.93180839162
s = 0.991019606076
Lua code had been updated, instead of starting s = a², it start at 0 (see post #7)
The reason is to produce as accurate E(m) - K(m) as possible.
lua> m = 0.25
lua> x, y = E.agm(1, sqrt(1-m))
lua> x, y
0.9318083916224482 -0.008980393923673079
lua> k = pi/(2*x)
lua> k, k*(y-m)/2 -- = K(m), E(m)-K(m)
1.685750354812596 -0.21828814547316888
RE: (28/48/50) Complete Elliptic Integrals - John Keith - 03-04-2024 09:15 PM
Thanks Albert, works fine now. Next step, comparison with my current programs for difficult values with |k| close to 1.
RE: (28/48/50) Complete Elliptic Integrals - Albert Chan - 03-05-2024 03:42 PM
(03-04-2024 06:39 PM)Albert Chan Wrote: Lua code had been updated, instead of starting s = a², it start at 0 (see post #7)
The reason is to produce as accurate E(m) - K(m) as possible.
I just realized I never posted lua update. Here it is
Code:
function E.agm(a, b)
local s, t = 0, 1
repeat
local b0, g = b, (a-b)/2
a, b = b+g, sqrt(a*b)
t = t + t
s = s - t*g*g
until b == b0 or not finite(b)
return b, s
endCode:
function E.K(m,c) -- K(m), E(m)-K(m), E(m)
if not c or m+c~=1 then c=1-m end -- c only a hint
if c==0 then return inf, -inf, 1 end
local x, y = E.agm(1, sqrt(c))
x = pi/(4*x)
return x+x, x*(y-m), x+x*(y+c)
endlua> eps = 1e-8
lua> E.K(1-eps, eps)
10.596634757087662 -9.596634706604487 1.000000050483174
Compare this to asymptote result for K(1-ε), E(1-ε)
lua> k0 = log(16/eps)/2
lua> k0 + eps/4*(k0-1) -- ≈ K(1-eps)
10.59663475708766
lua> 1 + eps/4*(2*k0-1) -- ≈ E(1-eps)
1.0000000504831736
RE: (28/48/50) Complete Elliptic Integrals - John Keith - 03-05-2024 04:05 PM
Interesting, seems a bit different than your Python program. I will try an RPL version and compare them.
Also, I found a bug in my translation of your Python program which would occasionally cause the loop to never terminate. I posted the corrected version above.
RE: (28/48/50) Complete Elliptic Integrals - John Keith - 03-09-2024 08:22 PM
Post #1 now contains new and updated programs. Kk and EKk have been improved based on Albert Chan's suggestion in post #7, resulting in cleaner code and improved accuracy for very small values of k. Ek has been removed due to redundancy; for Ek, use EKk DROP.
Two new programs have been added to compute the complete elliptic integral of the third kind and the elliptic nome q(x). See post #1 for explanation and links.
RE: (28/48/50) Complete Elliptic Integrals - John Keith - 03-14-2024 07:16 PM
Here are four short programs to compute the derivatives of the elliptic integrals at k. For references see the Wikipedia pages here and here.
The programs dKk, dEk and dNOME require k on the stack. dΠdk requires the numbers n on level 2 and k on level 1. I am listing the programs as a directory here, and they need to be in the path of the programs in post #1, several of which are called by these programs.
Code:
DIR
dKk
\<< 1 OVER SQ - OVER EKk 3 PICK * - ROT ROT * /
\>>
dEk
\<< DUP EKk - SWAP /
\>>
d\PIdk
\<<
IF OVER ABS
THEN DUP2 \PInk ROT ROT DUP EKk DROP
OVER SQ 1 - /
4 ROLL + OVER *
ROT ROT SQ - /
ELSE SWAP DROP dKk
END
\>>
dNOME
\<< DUP ENOME \pi SQ \->NUM *
SWAP 1 OVER SQ - OVER Kk SQ
ROT 2 * * * /
\>>
ENDAnd again, the same programs optimized for the HP 49 and 50.
Code:
DIR
dKk
\<< 1. OVER SQ - OVER EKk PICK3 * - UNROT * /
\>>
dEk
\<< DUP EKk - SWAP /
\>>
d\PIdk
\<<
IF OVER ABS
THEN DUP2 \PInk UNROT DUP EKk DROP
OVER SQ 1. - /
4. ROLL + OVER *
UNROT SQ - /
ELSE NIP dKk
END
\>>
dNOME
\<< DUP ENOME \pi SQ \->NUM *
SWAP 1. OVER SQ - OVER Kk SQ
ROT 2. * * * /
\>>
END