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

[Thomas Okken]

(06-19-2020 01:01 AM)DM48 Wrote:  In Werner's program below, is there anyway to make LBL "FX" into a local label? I keep trying to turn FX into 11 but this gives an error.

No, that is not possible. The PGMINT function requires a global label as its argument.

[DM48]

(06-19-2020 01:46 AM)Thomas Okken Wrote:  

(06-19-2020 01:01 AM)DM48 Wrote:  In Werner's program below, is there anyway to make LBL "FX" into a local label? I keep trying to turn FX into 11 but this gives an error.

No, that is not possible. The PGMINT function requires a global label as its argument.

Thanks Thomas. I was afraid that would be the case.

[Werner]

No need to be afraid..

Code:

00 { 96-Byte Prgm }
01▸LBL 10
02 VARMENU "BORE"
03 STOP
04 PGMINT "BORE"
05 1ᴇ-5
06 STO "ACC"
07 CLX
08 STO "LLIM"
09 RCL "C"
10 STO "ULIM"
11 INTEG "X"
12 4
13 ×
14▸LBL "BORE"
15 MVAR "A"
16 MVAR "B"
17 MVAR "C"
18 FC? 46 @ not integrating?
19 GTO 10
20 RCL "B"
21 XEQ 01
22 RCL "A"
23 XEQ 01
24 -
25 RCL "C"
26 XEQ 01
27 ×
28 RTN
29▸LBL 01
30 X↑2
31 RCL "X"
32 X↑2
33 -
34 SQRT
35 END

Cheers, Werner

[DM48]

Werner! You never cease to blow my mind. I assumed, incorrectly I see, that the first Global Label was the name of the program. Thank you for showing me this doesn't have to be the case.

[Werner]

There’s only one global label here, but it doesn’t have to be the first line in the program; it saved three bytes this way! Old habits die hard ;-)
Werner

[Thomas Okken]

Strictly speaking, programs don't have names. They can contain zero or more global labels, but most functions that take labels as arguments refer to that specific label, not to the program as a whole. This applies to GTO, XEQ, PGMINT, VARMENU, and most others.

The only exceptions, where a function does refer to an entire program, meaning, everything following the previous program's END (or the start of memory, in the case of the first program in memory) until the END, are the CLP and PRP functions, and for those functions, it does not matter which label within a program you choose.

[DM48]

(06-19-2020 10:42 PM)Werner Wrote:  There’s only one global label here, but it doesn’t have to be the first line in the program; it saved three bytes this way! Old habits die hard ;-)
Werner

I find it very slick you used the VARMENU with MVAR. Learning something new everyday.

[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

[DM48]

24 decimal precision (did I count correctly)? Wow!

[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

[Werner]

(06-20-2020 04:29 PM)DM48 Wrote:  24 decimal precision (did I count correctly)? Wow!

No, 34 digits!
Werner

[Albert Chan]

(06-06-2020 08:33 PM)Albert Chan Wrote:  

(06-02-2020 08:31 AM)Werner Wrote:  Question: wouldn't the formula need to be

bore hole volume = hv(b,c) - hv(max(a,c),min(a,c))

so that it would work for a<c as well?

Amazingly, the answer is NO !
...
The real part returns the same volume !

We can prove this with Legendre's relation, 19.7.3
Note: the prove use elliptical parameter m instead of modulus, k, m = k²

Assume d < D = 1, so that m = (d/D)² = d² < 1
(we can scale it back later, by factor D³)

3*HV(d,1) = m*(E+K) + (E-K) = (m+1)*E + (m-1)*K

With order flipped, using notation K(1/m) = K(m) , K(1-m) = K'(m)

3*HV(1,d)
= d³ * 3*HV(1/d,1)
= d³ * ((1/m+1)*E + (1/m-1)*K)
= (m+1) * (d*E) + (m-m²) * (K/d)
= (m+1) * (E - (1-m)*K ± 1j * (E' - m*K')) + (m-m²) * (K ∓ 1j * K')

Re(3*HV(1,d)) = (m+1)*E + (m²-1)*K + (m-m²)*K = 3*HV(d,1)

Adding the special case, HV(1,1) = 2/3, and removed assumptions d < D = 1:

hole volume = Re(HV(d,D)) = Re(HV(D,d))