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

[Albert Chan]

Hi, Werner

I added your suggested changes. Thanks.

(06-03-2020 06:45 PM)Albert Chan Wrote:  There is still an issue of AGM convergence.
For complex numbers, AGM might not converge *exactly*.
For example, on Free42, 2 XEQ "EK" crash with "Out-of-Range"

For EK(2), this involved calculating:

b = sqrt(1-m) = 1j
x1, y1 = agm2(1+b, 2*sqrt(b)) = agm2(1+1j, sqrt(2)*(1+1j))

We can scale away 1+1j, and get x1,y1, like this:

x2, y2 = agm2(1, sqrt(2))             ≈ (1.1981, 0.91389)
x1 = (1+1j) * x2                          ≈ 1.1981 + 1.1981j
y1 = (1+1j)**2 * y2 = 2j * y2      ≈ 1.8278j

→ K(2) = pi / x1                           ≈ 1.3110 - 1.3110j
→ E(2) = K(2) * y1/4                    ≈ 0.59907 + 0.59907j

We have y1 as purely imaginary complex number.
However, with rounding errors, directly go for y1 might have tiny real part.
Since we test for ACC equality, going directly for y1 might not work.

This is my idea ... by going up the agm2 chain.

x0, y0 = agm2(1,b),       b ≠ 0

x0 = agm(1,b) = agm((1+b)/2, sqrt(b)) = x1 / 2

y1 = (1+b)^2                - 2*((1+b)/2 - sqrt(b))^2 - 4*gap2^2 - 8*gap3^2 - ...
2 y0 = 2*1^2 - (1-b)^2 - 2*((1+b)/2 - sqrt(b))^2 - 4*gap2^2 - 8*gap3^2 - ...

→ y1 - 2 y0 = (1+b)^2 - 2 + (1-b)^2 = 2 b^2
→ y1 = 2*(y0 + b^2)

x0, y0 = agm2(1, 1j)                    ≈ (0.59907 + 0.59907j, 1.0000 + 0.91389j)

→ K(2) = pi / (2*x0)                    ≈ 1.3110 - 1.3110j
→ E(2) = K(2) * (y0 + b^2)/2      ≈ 0.59907 + 0.59907j

This is my new version, getting EK(m) using AGM(1,b), b = sqrt(1-m)

Code:

00 { 181-Byte Prgm }
01▸LBL "BORE"
02 LSTO "c"
03 R↓
04 LSTO "b"
05 R↓
06 XEQ 03
07 LSTO "t"
08 RCL "b"
09 XEQ 03
10 RCL "t"
11 RTN
12▸LBL 03
13 RCL "c"
14 X>Y?
15 X<>Y
16 X↑2
17 LSTO "d"
18 X<>Y
19 LSTO "D"
20 X↑2
21 ÷
22 XEQ "EK"
23 RCL- ST Y
24 X<>Y
25 LASTX
26 +
27 RCL× "d"
28 X<>Y
29 RCL "D"
30 X↑2
31 ×
32 +
33 RCL× "D"
34 3
35 ÷
36 RTN
37▸LBL "EK"
38 1
39 -
40 ABS
41 X=0?
42 GTO 00
43 1
44 LASTX
45 +/-
46 LSTO "k"
47 SQRT
48 XEQ "AGM"
49 PI
50 X<>Y
51 STO+ ST X
52 ÷
53 X<>Y
54 RCL+ "k"
55 RCL× ST Y
56 2
57 ÷
58 RTN
59▸LBL 00
60 99
61 1
62 RTN
63▸LBL "AGM"
64 LSTO "B"
65 ABS
66 SIGN
67 +/-
68 LSTO "T"
69 R↓
70 LSTO "A"
71 X↑2
72▸LBL 02
73 RCL "A"
74 RCL× "B"
75 SQRT
76 X<> "B"
77 RCL- "A"
78 2
79 STO× "T"
80 ÷
81 STO+ "A"
82 X↑2
83 RCL× "T"
84 RCL+ ST Y
85 X≠Y?
86 GTO 02
87 RCL "B"
88 END

2 XEQ "EK"     → E ≈ 0.59907 + 0.59907i,     K ≈ 1.3110 - 1.3110i

Update: replace EK(m = (d/D)^2) as EK(m = d^2/D^2)
Possibly more accurate m, and saved 1 byte of code.