(28/48/50) Complete Elliptic Integrals

[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