(28/48/50) Complete Elliptic Integrals
02-26-2024, 08:14 PM (This post was last modified: 03-11-2024 08:48 PM by John Keith.)
Post: #1
 John Keith Senior Member Posts: 1,028 Joined: Dec 2013
(28/48/50) Complete Elliptic Integrals
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   \>> END

And 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   \>> END
02-27-2024, 12:08 AM (This post was last modified: 02-27-2024 12:23 AM by Thomas Klemm.)
Post: #2
 Thomas Klemm Senior Member Posts: 2,100 Joined: Dec 2013
RE: (28/48/50) Complete Elliptic Integrals
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
02-27-2024, 08:53 PM
Post: #3
 John Keith Senior Member Posts: 1,028 Joined: Dec 2013
RE: (28/48/50) Complete Elliptic Integrals
(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.
02-28-2024, 11:21 AM
Post: #4
 Gjermund Skailand Member Posts: 78 Joined: Dec 2013
RE: (28/48/50) Complete Elliptic Integrals
There is also my old post in the old forum which you may find interesting

br Gjermund
02-28-2024, 03:57 PM
Post: #5
 floppy Senior Member Posts: 363 Joined: Apr 2021
RE: (28/48/50) Complete Elliptic Integrals
Using AGM and MAGM is perhaps quicker?
https://stackoverflow.com/questions/4231...0#69606140

HP71 4TH/ASM & Multimod, HP41CV/X & Nov64d, PILBOX, HP-IL 821.62A & 64A & 66A, Deb11 64b-PC & PI2 3 4 w/ ILPER, VIDEO80, V41 & EMU71, DM41X, HP75D
02-28-2024, 09:54 PM
Post: #6
 John Keith Senior Member Posts: 1,028 Joined: Dec 2013
RE: (28/48/50) Complete Elliptic Integrals
(02-28-2024 03:57 PM)floppy Wrote:  Using AGM and MAGM is perhaps quicker?
https://stackoverflow.com/questions/4231...0#69606140

For the perimeter of the ellipse perhaps, but these programs compute the elliptic integrals which may be used for other purposes.
02-29-2024, 04:31 AM
Post: #7
 Albert Chan Senior Member Posts: 2,677 Joined: Jul 2018
RE: (28/48/50) Complete Elliptic Integrals
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
02-29-2024, 09:23 PM
Post: #8
 John Keith Senior Member Posts: 1,028 Joined: Dec 2013
RE: (28/48/50) Complete Elliptic Integrals
Thanks, Albert. I will try your ideas and Werner's, although it will complicate the program- RPL has neither GOTO nor BREAK.
02-29-2024, 11:43 PM
Post: #9
 Thomas Klemm Senior Member Posts: 2,100 Joined: Dec 2013
RE: (28/48/50) Complete Elliptic Integrals
(02-29-2024 09:23 PM)John Keith Wrote:  RPL has neither GOTO nor BREAK.

Cf. RPL equivalents to BREAK, CYCLE, EXIT, etc.?
03-04-2024, 05:37 PM
Post: #10
 John Keith Senior Member Posts: 1,028 Joined: Dec 2013
RE: (28/48/50) Complete Elliptic Integrals
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.
03-04-2024, 06:11 PM (This post was last modified: 03-05-2024 03:53 PM by John Keith.)
Post: #11
 John Keith Senior Member Posts: 1,028 Joined: Dec 2013
RE: (28/48/50) Complete Elliptic Integrals
(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.
03-04-2024, 06:39 PM
Post: #12
 Albert Chan Senior Member Posts: 2,677 Joined: Jul 2018
RE: (28/48/50) Complete Elliptic Integrals
(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
03-04-2024, 09:15 PM
Post: #13
 John Keith Senior Member Posts: 1,028 Joined: Dec 2013
RE: (28/48/50) Complete Elliptic Integrals
Thanks Albert, works fine now. Next step, comparison with my current programs for difficult values with |k| close to 1.
03-05-2024, 03:42 PM
Post: #14
 Albert Chan Senior Member Posts: 2,677 Joined: Jul 2018
RE: (28/48/50) Complete Elliptic Integrals
(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 end
Code:
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) end

lua> 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
03-05-2024, 04:05 PM
Post: #15
 John Keith Senior Member Posts: 1,028 Joined: Dec 2013
RE: (28/48/50) Complete Elliptic Integrals
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.
03-09-2024, 08:22 PM
Post: #16
 John Keith Senior Member Posts: 1,028 Joined: Dec 2013
RE: (28/48/50) Complete Elliptic Integrals
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.
03-14-2024, 07:16 PM
Post: #17
 John Keith Senior Member Posts: 1,028 Joined: Dec 2013
RE: (28/48/50) Complete Elliptic Integrals
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 * * * /   \>> END

And 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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)