(28/48/50) Complete Elliptic Integrals

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