HP42s first major program (Double Integral) Best way to approach?

[Albert Chan]

(06-13-2020 10:16 PM)Albert Chan Wrote:  Inspired by Werner's AGM code thread, I fixed my AGM2 convergence issue.

It now test convergence of sucessive GM's, eliminated convergence issue of AGM(x,0).

To improve agm2 accuracy, I redefined agm2 returns:

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

Below is updated EK based HV, and K/AGM code

Code:

00 { 65-Byte Prgm }
01▸LBL "HV"
02 LSTO "c"
03 X<>Y
04 LSTO "b"
05 R↓
06 XEQ 02       ; HV(a,c)
08 RCL "c"      ; HV(b,c)
12▸LBL 02
13 X>Y?
14 X<>Y
15 X↑2
16 X<>Y
17 LSTO "D"
18 X↑2
19 STO× "D"     ; D := D^3
20 ÷
21 XEQ 03       ; HV1(m)
22 RCL× "D"
23 ENTER
24 X<> "b"
23 RTN
24▸LBL 03       ; HV1(m) = (m*(E+K) + (E-K))/3
25 SIGN         ; m = (d/D)^2 >= 0, sign = 1
26 LASTX
27 X=Y?
28 GTO 04       ; HV1(1) = 2/3
29 XEQ "K"      ; K    E-K  m   m
30 RCL+ ST Y
31 LASTX
32 +            ; E+K  E-K  m   m
33 R↑
34 ×
35 LBL 04
36 +
37 3
38 ÷
39 END

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

(06-08-2020 12:11 AM)Albert Chan Wrote:  HV(1,1000) = 785.3980652226656130844841050849260, error = -2 ULP
BORE       = 785.3980652226656130844841051786603, error = -937345 ULP
Mathemtaica: 785.3980652226656130844841050849257781988463583655713086050...

This updated version matched accuracy of my EKmc based HV.
With EK based HV(1,1000), I get:

785.3980652226656130844841050849257