(28/48/50) Complete Elliptic Integrals

+- 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 03: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
  \>>
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


RE: (28/48/50) Complete Elliptic Integrals - John Keith - 03-14-2024 02: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 * * * /
  \>>
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