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

[Albert Chan]

Hi, Werner

I made some improvements to BORE program ...

Code:

00 { 176-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 LSTO "d"
17 X<>Y
18 LSTO "D"
19 ÷
20 X↑2
21 XEQ "EK"
22 RCL- ST Y
23 X<>Y
24 LASTX
25 +
26 RCL "d"
27 X↑2
28 ×
29 X<>Y
30 RCL "D"
31 X↑2
32 ×
33 +
34 RCL× "D"
35 3
36 ÷
37 RTN
38▸LBL "EK"
39 1
40 X=Y?
41 GTO 00
42 -
43 LASTX
44 X<>Y
45 +/-
46 SQRT
47 STO+ ST Y
48 SQRT
49 STO+ ST X
50 XEQ "AGM"
51 PI
52 X<>Y
53 ÷
54 STO× ST Y
55 X<>Y
56 4
57 ÷
58 RTN
59▸LBL 00
60 99
61 X<>Y
62 RTN
63▸LBL "AGM"
64 LSTO "B"
65 RCL÷ "B"
66 +/-
67 LSTO "T"
68 R↓
69 LSTO "A"
70 X↑2
71▸LBL 02
72 RCL "A"
73 RCL× "B"
74 SQRT
75 X<> "B"
76 RCL- "A"
77 2
78 STO× "T"
79 ÷
80 STO+ "A"
81 X↑2
82 RCL× "T"
83 RCL+ ST Y
84 X≠Y?
85 GTO 02
86 RCL "B"
87 END

First, I followed Mathematica EllipticE, EllipticK, using parameter m instead of modulus k, m = k^2

With this change, for ellipse_perimeter(10,50), a = 50, e^2 = 1-0.2^2 = 0.96

0.96 XEQ "EK" 200 ×       → 210.100445397

We can do the same calculation with a = 10, e^2 = 1-5^2 = -24
(previous version can do this too, but with imaginary e)

-24 XEQ "EK" 40 ×          → 210.100445397

I noticed K(m) is relatively slow growing, when m approach 1.
Example: K(m = 1 - 1e-85) ≈ 99

Instead of adding the rule HV(c,c) = 2/3*c^3, I do E(1)=1, K(1)=99 (see LBL 00)
This kill 2 birds with 1 stone

I also removed the menu, with "BORE" getting a,b,c, directly from the stack.
To avoid catastrophic cancellation problems without knowing it, I leave HV values on the stack.

48 Enter 58 Enter 48
XEQ "BORE"             → Y=94954.7588409, X=73728
−                            → I=21226.7588409

This allowed removal of HV label: Example, for HV(29,12)

29 Enter 12 Enter
XEQ "BORE"             → 3208.03490414

I also removed the assumption SIGN giving 1, and add 1 explicitly.
This allowed the code to run correctly with complex numbers.

3 XEQ "EK"             → E ≈ 0.47522 + 1.0130j, K ≈ 1.0011 - 1.1714j
-1 SQRT XEQ "EK"   → E ≈ 1.6324 - 0.36922j, K ≈ 1.4213 + 0.29538j

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"