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

[Albert Chan]

(06-20-2020 04:23 PM)Albert Chan Wrote:  x, y = agm2(1, sqrt(1-m))
→ K = pi / (2*x)
→ E = K * (y + (1-m) + 1²) / 2 = K + K*(y-m)/2

We can also move up the agm2 chain, and do this (note, m=d²)

x0, y0 = agm2(1+d, 1-d)
→ K = pi / (2*x0)
→ E = K + K*y0/4

→ HV1 = (m*(E+K) + (E-K))/3 = pi/(24*x0) * (y0*(m+1) + 8*m)

We can get HV without E, K, not even sqrt, only agm2.
For example, HV(1, 1000)

lua> d, D = 1, 1000
lua> d2, D2 = d*d, D*D
lua> x, y = agm2(D+d, abs(D-d))
lua> pi/(24*x) * (y*(d2+D2) + 8*d2*D2)       -- HV(d,D)
785.3980652226655

Note the symmetry of HV formula, it does not require d ≤ D (no sorting )

Code:

00 { 71-Byte Prgm }
01▸LBL "HV"
02 LSTO "c"
03 R↓
04 LSTO "b"
05 R↓
06 XEQ 03       ; HV(a,c)
07▸LBL 03       ; HV = pi/x/3 * (y*(dd+DD)/8 + dd*DD)
08 LSTO "d"
09 STO ST Y
10 RCL "c"
11 STO+ ST Z    ; D     d    D+d
12 X=Y?
13 GTO 04       ; HV(d,d) = 2/3*d^3
14 -
15 ABS
16 XEQ "AGM"    ; agm(D+d,|D-d|) -> (x,y)
17 X<> "d"
18 X↑2
19 STO ST Z     ; dd    y   dd
20 RCL "c"
21 X↑2
22 STO× ST T    ; DD    dd  y   dd*DD
23 +
24 ×
25 8
26 ÷
27 +
28 PI
29 RCL "d"
30 1/X
31▸LBL 04
32 ×
33 ×
34 3
35 ÷
36 ENTER
37 X<> "b"
38 END