(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