(28/48/50) Complete Elliptic Integrals

[John Keith]

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

[Thomas Klemm]

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

[John Keith]

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

[Gjermund Skailand]

There is also my old post in the old forum which you may find interesting
: http://www.hpmuseum.org/cgi-sys/cgiwrap/...ead=245627

br Gjermund

[floppy]

Using AGM and MAGM is perhaps quicker?
https://stackoverflow.com/questions/4231...0#69606140

[John Keith]

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

[Albert Chan]

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 = -Σ(2^k (½gap_k)² , 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

[John Keith]

Thanks, Albert. I will try your ideas and Werner's, although it will complicate the program- RPL has neither GOTO nor BREAK.

[Thomas Klemm]

(02-29-2024 09:23 PM)John Keith Wrote:  RPL has neither GOTO nor BREAK.

Cf. RPL equivalents to BREAK, CYCLE, EXIT, etc.?

[John Keith]

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.

[John Keith]

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

[Albert Chan]

(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

[John Keith]

Thanks Albert, works fine now. Next step, comparison with my current programs for difficult values with |k| close to 1.

[Albert Chan]

(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

[John Keith]

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.

[John Keith]

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.

[John Keith]

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