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

[Werner]

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?

Werner

[Albert Chan]

(06-02-2020 08:23 AM)Werner Wrote:  Continuous improvement.. when I rewrite agm as ...

My version tried to avoid last sqrt(), but your version had last sqrt() evaluated.
You might as well return that instead of arithmetic mean.

function agm2(a, b)
local t = -1
local s1 = a*a
repeat
local s = s1
local k = (b-a)/2
t = t + t
b = sqrt(a*b)
a = a + k
s1 = s + t*k*k
until s == s1
return b, s1
end

(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?

Same reasoning goes for b<c

A simple fix is to sort the hv arguments, such that hv(a,c) = hv(c,a)

[Werner]

okay, final update then:
will now work when b<c and/or a<c, as well; and a few more optimisations.

Code:

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

Werner

[DM48]

When a = 48, b = 58, and c = 48


BORE (Integration)
Answer = (one side) 10,613 * 2 => (both sides) 21,226.76 units^3


BORE (Werners final program)
I get a Divde by 0 error on Line 59. Shouldn't this only return an error (non-real result) when c > a?

[Massimo Gnerucci]

48 is not RPN. Try again.

SCNR...

[Albert Chan]

Hi, DM48

HV(c,c) involve calculating K(1) = ∞, thus Divide by 0 error.

The fix required another rule: HV(c,c) = 2/3*c^3

For a=48, b=58, c=48:

HV(58,48) ≈ 94954.7588409
HV(48,48) = 2/3*48^3 = 73728

Difference ≈ 21226.7588409

[Werner]

Easy fix

Code:

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

Werner

[DM48]

(06-03-2020 05:31 AM)Werner Wrote:  Easy fix

Code:

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

Werner

Werner, testing the fix, it does work at c = 48.

Since I only need half the volume, I use to slide in a 2 / after line 43 in this program. That returns a large negative result. Do you have a suggestion where I should place 2 / so it performs the division on the final result?

[Werner]

Yes, between lines 13 and 14.
Werner

[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"

[Werner]

Just one more thing:
EK has to test for input 1 and (1,0), so it has to be (never mind the step numbers):

Code:

42▸LBL "EK"
43 1
44 -
45 ABS
46 X=0?
47 GTO 00
48 1
49 LASTX
50 +/-
51 SQRT
52 STO+ ST Y
53 SQRT
54 STO+ ST X
55 XEQ "AGM"
56 PI
57 X<>Y
58 ÷
59 STO× ST Y
60 X<>Y
61 4
62 ÷
63 RTN
64▸LBL 00
65 99
66 1
67 RTN
68▸LBL "AGM"
69 LSTO "B"
70 ABS
71 SIGN
72 +/-
73 LSTO "T"

I added my start of AGM as well, shorter, and works for B=0

Cheers, Werner

[DM48]

Are these lines meant to replace lines 38 - 67 of Albert's code?

(06-04-2020 02:00 PM)Werner Wrote:  Just one more thing:
EK has to test for input 1 and (1,0), so it has to be (never mind the step numbers):

Code:

42▸LBL "EK"
43 1
44 -
45 ABS
46 X=0?
47 GTO 00
48 1
49 LASTX
50 +/-
51 SQRT
52 STO+ ST Y
53 SQRT
54 STO+ ST X
55 XEQ "AGM"
56 PI
57 X<>Y
58 ÷
59 STO× ST Y
60 X<>Y
61 4
62 ÷
63 RTN
64▸LBL 00
65 99
66 1
67 RTN
68▸LBL "AGM"
69 LSTO "B"
70 ABS
71 SIGN
72 +/-
73 LSTO "T"

I added my start of AGM as well, shorter, and works for B=0

Cheers, Werner

[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.

[Albert Chan]

(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 !

Example: a, b, c = 48, 58, 68

The drill bit is bigger than the pipe, and will cut it off.

48 Enter 58 Enter 68 XEQ "BORE"

→ HV(c,b) = 161316.6882015429652660165715495419
→ HV(c,a) = 114811.9333705716751263443462068749
Difference =   46504.754830971290139672225342667

My latest code update should handle complex numbers.
Now, remove the HV arguments sorting steps 14, 15, and redo above example

→ HV(c,b) = 161316.688201542965266016571549542 + 4944.10953174134265265013475228068i
→ HV(c,a) = 114811.9333705716751263443462068748 + 18153.27669892468194546755520844514i
Difference =   46504.7548309712901396722253426672 - 13209.16716718333929281742045616446i

The real part returns the same volume !
Bonus: imaginery part warned the user that pipe is cut off

[Albert Chan]

This version uses my implementation of EKmc(m) = E(m) - c, K(m) - c, where c = Pi/2

I added HV(d,D) on top of EKmc(m) to test its accuracy.

HV(d, D) = D/3 (d^2 (E+K) + D^2 (E-K))
              = D/3 (d^2 ((E-c)+(K-c) + Pi) + D^2 ((E-c)-(K-c)))

Note: this is a recursive algorithm, requiring Free42 LSTO command.

Code:

00 { 167-Byte Prgm }
01▸LBL "HV"
02 LSTO "c"
03 X<>Y
04 LSTO "b"
05 R↓
06 XEQ 03       ; HV(c,a)
07 X<> "b"
08 RCL "c"
09 XEQ 03       ; HV(c,b)
10 RCL "b"
11 RTN
12▸LBL 03       ; HV(d,D)
13 X>Y?         ; = D/3*(dd*(e+k) + DD*(e-k))
14 X<>Y         ; = D/3*((dd+DD)*(e-c) + (dd-DD)*(k-c) + dd*(2c))
15 X↑2
16 LSTO "d"     ; d := d^2
17 X<>Y
18 LSTO "D"
19 X=Y?
20 GTO 01       ; HV(d,d) = 2/3*d^3
21 X↑2
22 ÷            ; m = (dd)/(DD)
23 XEQ "EKmc"   ; e-c    k-c
24 RCL "d"
25 RCL "D"
26 X↑2
27 -
28 STO× ST Z    ; dd-DD  e-c    (dd-DD)*(k-c)
29 X<> ST L
30 RCL+ "d"     ; dd+DD  e-c    (dd-DD)*(k-c)
31 ×
32 +
33 PI
34 RCL× "d"
35▸LBL 00
36 +
37 RCL× "D"
38 3
39 ÷
40 RTN
41▸LBL 01       ; (D) -> 2/3*D^3
42 X↑2
43 ENTER
44 GTO 00
45▸LBL "EKmc"   ; E(m)-c, K(m)-c, where c=PI/2
46 ENTER
47 ABS
48 1
49 +
50 LASTX        ; 1  1+|m|  m
51 X=Y?
52 GTO 02       ; EKmc base case
53 STO ST Y
54 RCL- ST Z    ; 1-m  1    m
55 SQRT
56 LSTO "b"     ; b    1   m
57 +
58 X↑2
59 ÷
60 LSTO "h"     ; h = (1-b)/(1+b) = m/(1+b)^2
61 X↑2
62 XEQ "EKmc"   ; EKmc(h*h) = e, k
63 RCL× "b"     ; EKmc(m) = (hk-k) + (e+be) - bhc, (hk+k) + hc
64 LASTX
65 +
66 X<>Y
67 RCL× "h"
68 LASTX
69 RCL+ ST Y
70 X<>Y
71 LASTX
72 -
73 R↑
74 +
75 PI
76 2
77 ÷
78 RCL× "h"
79 STO+ ST Z
80 RCL× "b"
81 -
82 RTN
83▸LBL 02
84 PI           ; PI    1   1+|m|   m
85 R↑
86 ×
87 8
88 ÷
89 ENTER
90 +/-          ; -pi/8*m   pi/8*m
91 END

Let's try HV(1,1000):

HV(1,1000) = 785.3980652226656130844841050849260, error = -2 ULP
BORE       = 785.3980652226656130844841051786603, error = -937345 ULP
Mathemtaica: 785.3980652226656130844841050849257781988463583655713086050...

* Mathematica: HV = Pi/4 d^2 D Hypergeometric2F1[-1/2,1/2,2,(d/D)^2] /. {d->1, D->1000}

Update: HV now has same functionality of BORE, and better accuracy.
Formula changed to avoid subtraction cancellation, when m = (d/D)² ≈ 1

HV(d, D) = D/3*((dd+DD)*(e-c) + (dd-DD)*(k-c) + dd*(2c))
→ HV1(m) = 1/3*((m+1)*(e-c) + (m-1)*(k-c) + m*pi)

First 2 terms both negative, with -sum less than 1/2 of the third.
We know this because 2/3 ≤ HV1/m ≤ pi/4, and (third term)/m = pi/3

→ (third term)/HV1 = [4/3, pi/2] ≈ [1.3333, 1.5708]
→ -(sum of first 2 terms)/(third term) ≈ [0.25000, 0.36338]

[Albert Chan]

Here is EKmc() non-recursive Python code.
Note: code preceded by "from math import sqrt, pi", for complex number support, use cmath.

Code:

def EKmc(m, verbal=False, c=pi/2):
    H = []          # non-recursive version, by collecting h's
    while abs(1+m)-1: h = m/(1+sqrt(1-m))**2; H.append(h); m = h*h
    k = 0.25*c*m
    e = -k
    for h in reversed(H):
        m = h*h
        if verbal: print 'm =',m,'\n\tE(m)-c =',e,'\n\tK(m)-c =',k
        e = (e+e-k + m*k + (m-h)*c) / (1+h)        
        k = k + h*k + h*c
    return e, k

Try HV(1,1000) again. With d=1, we can simplify the HV formula.

>>> D = 1000
>>> e, k = EKmc(1/(D*D), verbal=True)
m = 2.44141113282e-28
      E(m)-c = -9.58739909907e-29
      K(m)-c = 9.58739909907e-29
m = 6.25000625001e-14
      E(m)-c = -2.45437171499e-14
      K(m)-c = 2.45437171499e-14
>>> e, k      # m = 1e-6
(-3.9269915532983253e-07, 3.9269930259211089e-07)

>>> D/3 * ((e+k) + D*D*(e-k) + pi)       # = HV(1,D)
785.39806522266554

EKmc has the speed of AGM2 method, but avoided catastrophic cancellation issue.
It also avoided AGM2 (complex number) convergence issue.

With K(m), we can recover AGM's: AGM(1, b) * K(1-b*b) = pi/2

[DM48]

Impressive Albert. I'll be chewing through this for a while. I need to get my programming skills up to speed.

[Albert Chan]

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).
Test for EK(m) for abs(m-1)=0 is now shifted to test returns of AGM2.

To simplify further , I make this version dimensionless, m = (d/D)^2

→ HV(d,D) = HV1(m) * D³

Code:

00 { 128-Byte Prgm }
01▸LBL "HV1"    ; HV1(m) = (m(e+k) + (e-k))/3
02 LSTO "m"     ; m = (d/D)^2
03 XEQ "EK"
04 RCL- ST Y
05 X<>Y
06 LASTX
07 +
08 RCL× "m"
09 +
10 3
11 ÷
12 RTN
13 ▸LBL "EK"    ; (m) -> E(m), K(m)
14 1
15 -
16 LASTX
17 X<>Y
18 +/-          ; 1-m   1
19 SQRT
20 STO+ ST Y
21 SQRT
22 STO+ ST X
23 XEQ "AGM"    ; x     y
24 ABS
25 X=0?
26 GTO 00
27 PI
28 LASTX        ; x     pi    abs(x)  y
29 ÷
30 STO× ST T
31 R↑
32 4
33 ÷
34 RTN
35▸LBL 00       ; K(1) assumed 99
36 99
37 1
38 RTN
39▸LBL "AGM"    ; (a,b) -> (agm, acc)
40 X<>Y
41 LSTO "A"
42 X↑2
43 LSTO "S"
44 1
45 LSTO "T"
46 X<> ST Z
47▸LBL 02       ; b
48 ENTER
49 RCL× "A"     ; ab    b
50 LASTX
51 RCL- "A"     ; b-a   ab    b
52 2
53 STO× "T"
54 ÷            ; k     ab    b     b
55 STO+ "A"
56 X↑2
57 RCL× "T"
58 STO- "S"
59 R↓
60 SQRT         ; GM    b     b     tkk
61 X≠Y?
62 GTO 02
63 RCL "S"
64 X<>Y         ; GM    S     b     b
65 END

[Albert Chan]

(05-27-2020 08:27 PM)DM48 Wrote:  I am looking to make this the fastest calculation possible with at least five, preferably six decimals.

Previous post HV1(m), if you plot it (m = 0 to 1), it look like a straight line, HV1 ≈ 2/3*m

Using a trick I learned earlier, for asinc(x), using Mathematica:

In[1]:= <<Calculus`Pade`
In[2]:= hv1[m_] := ((m+1) EllipticE[m] + (m-1) EllipticK[m]) / 3
In[3]:= EconomizedRationalApproximation[hv1[m]/m, {m,{0,0.75},2,2}]

Results coded in Lua, for m ≤ 0.7, this get about 6 digits accuracy.

Code:

function hv1(m)
    local num = (m-5.375879977642193) * (m-1.7259819956617186)
    local den = (m-12.968817144843689) * (m-1.7752730551508353)
    return 1.948810875414419 * m * num / den
end

Adding back dimensions to test the result:

lua> function hv(d,D) return D^3 * hv1((d/D)^2) end
lua> v2 = hv(24, 58)
lua> v1 = hv(24, 48)
lua> print(v2, v1, v2-v1)
25664.275659138966       21013.02544172724       4651.250217411725

For reference, this is EKmc HV version

Free42: 48 Enter 58 Enter 24 XEQ "HV" -

   25664.27923311925502790987133042097
− 21013.03279650214793028432717211253
=   4651.24643661710709762554415830844

[DM48]

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.

My goal is not to see FX in the list of programs, make BORE self contained if you will. I assume this has to do with the XEQ or PGMINT call?

(05-28-2020 07:28 AM)Werner Wrote:  Or you could use this, and use it interactively:

Code:

00 { 98-Byte Prgm }
01▸LBL "BORE"
02 MVAR "A"
03 MVAR "B"
04 MVAR "C"
05 1ᴇ-5
06 STO "ACC"
07 CLX
08 STO "LLIM"
09 PGMINT "FX"
10▸LBL 10
11 VARMENU "BORE"
12 STOP
13 RCL "C"
14 STO "ULIM"
15 INTEG "X"
16 4
17 ×
18 GTO 10
19▸LBL "FX"
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