integral competition HP50g vs. DM42

[Thomas Okken]

(08-24-2020 10:30 AM)peacecalc Wrote:  The algorithmus can be a kind of gaussian quadrature, because that one doesn't use the limits on the left and right side.

Romberg and other method use the limits of the integration aera.

Free42 uses are modified Romberg method that doesn't use the endpoints. The code was originally written by Hugh Steers; this comment is from core_math1.cc:

Code:


/* approximate integral of `f' between `a' and `b' subject to a given
 * error. Use Romberg method with refinement substitution, x = (3u-u^3)/2
 * which prevents endpoint evaluation and causes non-uniform sampling.
 */

[Joe Horn]

(08-24-2020 04:28 PM)Thomas Okken Wrote:  Free42 uses are modified Romberg method that doesn't use the endpoints.

That sounds familiar: http://holyjoe.net/HP71/integral.htm

[peacecalc]

Hello all,

Thank you: @Joe Horn, @John Keith and @Thomas Okken!

It is known, why the algorithmus was not switched through the generations of calculators?

It is a kind of "never change a running system"?

Most interesting was for me the transformation from x to u avoiding misleading results for periodic functions.

[JMBaillard]

Hi !
here is an HP41 version that employs an ascending series for small x
and a complex continued fraction otherwise:

01 LBL "CSX"
02 STO 01
03 ABS
04 PI
05 SQRT
06 X>Y?
07 GTO 04
08 DEG
09 *
10 24
11 STO 02
12 CLX
13 2
14 /
15 STO 03
16 ENTER
17 LBL 01
18 X<>Y
19 CHS
20 STO Z
21 X^2
22 RCL Y
23 X^2
24 +
25 DSE 02
26 X<0?
27 GTO 02
28 RCL 02
29 X<>Y
30 ST+ X
31 /
32 ST* Z
33 *
34 RCL 03
35 ST+ Z
36 +
37 GTO 01
38 LBL 02
39 PI
40 SQRT
41 *
42 ST/ Z
43 /
44 R-P
45 RCL 01
46 X^2
47 90
48 *
49 R^
50 -
51 X<>Y
52 CHS
53 P-R
54 1
55 +
56 STO Z
57 X<>Y
58 ST+ Z
59 -
60 2
61 ST/ Z
62 /
63 GTO 05
64 LBL 03
65 X<> 00
66 PI
67 2
68 ST+ 02
69 /
70 RCL 01
71 X^2
72 *
73 X^2
74 *
75 RCL 02
76 ENTER
77 DSE X
78 *
79 /
80 CHS
81 STO 00
82 RCL 02
83 ST+ X
84 ISG X
85 CLX
86 /
87 X<>Y
88 ST+ Y
89 X#Y?
90 GTO 03
91 RTN
92 LBL 04
93 CLX
94 STO 02
95 X<>Y
96 STO 00
97 ENTER
98 XEQ 03
99 STO 03
100 1
101 STO 02
102 RCL 01
103 ABS
104 3
105 Y^X
106 PI
107 *
108 2
109 /
110 STO 00
111 3
112 /
113 ENTER
114 XEQ 03
115 RCL 03
116 LBL 05
117 RCL 01
118 SIGN
119 ST* Y
120 ST* Z
121 RDN
122 END

Example:

1.2 XEQ "CSX" >>>> C(x) = 0.715437723
X<>Y S(x) = 0.623400918

3.9 R/S >>>> C(x) = 0.422332710
X<>Y S(x) = 0.475202402

Best regards.

[pinkman]

On my Prime G1 rev C I defined C(X) and S(X) with the Define key, and created two small programs looping on C(3.9) or S(3.9) and counting the TICKS difference.

The results:
Average C(3.9) time: 74.05ms
Average S(3.9) time: 73.67ms

[peacecalc]

Hello JMBaillard,

I converted your program "CSX" so, that's working for the DM42 (maybe for the HP42, too, but I've not this calc).

Code:


00 { 166-Byte Prgm }
01▸LBL "CSX"
02 STO 01
03 ABS
04 PI
05 SQRT
06 X>Y?
07 GTO 04
08 DEG
09 ×
10 24
11 STO 02
12 CLX
13 2
14 ÷
15 STO 03
16 ENTER
17▸LBL 01
18 X<>Y
19 +/-
20 STO ST Z
21 X↑2
22 RCL ST Y
23 X↑2
24 +
25 DSE 02
26 X<0?
27 GTO 02
28 RCL 02
29 X<>Y
30 STO+ ST X
31 ÷
32 STO× ST Z
33 ×
34 RCL 03
35 STO+ ST Z
36 +
37 GTO 01
38▸LBL 02
39 PI
40 SQRT
41 ×
42 STO÷ ST Z
43 ÷
44 →POL
45 RCL 01
46 X↑2
47 90
48 ×
49 R↑
50 -
51 X<>Y
52 +/-
53 →REC
54 1
55 +
56 STO ST Z
57 X<>Y
58 STO+ ST Z
59 -
60 2
61 STO÷ ST Z
62 ÷
63 GTO 05
64▸LBL 03
65 X<> 00
66 PI
67 2
68 STO+ 02
69 ÷
70 RCL 01
71 X↑2
72 ×
73 X↑2
74 ×
75 RCL 02
76 ENTER
77 DSE ST X
78 ×
79 ÷
80 +/-
81 STO 00
82 RCL 02
83 STO+ ST X
84 ISG ST X
85 CLX
86 ÷
87 X<>Y
88 STO+ ST Y
89 X≠Y?
90 GTO 03
91 RTN
92▸LBL 04
93 CLX
94 STO 02
95 X<>Y
96 STO 00
97 ENTER
98 XEQ 03
99 STO 03
100 1
101 STO 02
102 RCL 01
103 ABS
104 3
105 Y↑X
106 PI
107 ×
108 2
109 ÷
110 STO 00
111 3
112 ÷
113 ENTER
114 XEQ 03
115 RCL 03
116▸LBL 05
117 RCL 01
118 SIGN
119 STO× ST Y
120 STO× ST Z
121 R↓
122 END

Some commands are different: f. i. STO+ X is replaced by STO+ ST X;
Have fun!

[grsbanks]

(08-30-2020 10:51 AM)peacecalc Wrote:  I converted your program "CSX" so, that's working for the DM42 (maybe for the HP42, too, but I've not this calc).

Actually, the encoder over on the SwissMicros website accepts both syntaxes but converts to HP-42S syntax. Quoting from https://www.swissmicros.com/dm42/decoder/#tabs-2

Quote:Many of the characters used in HP-42S/DM42 programs cannot be entered directly with a keyboard. In order to help with this, we provide several special "codes" that you can use and there are also easy-to-type mnemonics for instructions that use them. This encoder also accepts HP-41 style instructions, eg. you can type "ST/ 10" instead of "STO÷ 10" or "RCL T" instead of "RCL ST T".

[JMBaillard]

Hi !
The HP41 works with 10 digits,
but to get 12-digit precision,
it's preferable to replace line 10 by
50 ( instead of 24 ).
Best regards.

[John Keith]

(08-23-2020 07:38 PM)Albert Chan Wrote:  >>> k = 2/pi**0.5
>>> erf = lambda x: k*x * hyp1f1(1/2, 3/2, -x*x)
>>> CiS = lambda x, erf=erf: (1+1j)/2 * erf((1-1j)/k * x)
>>> CiS(3.9)
(0.42233270272270274 + 0.47520241152278453j)

CiS(3.9) called erf(z = (1-1j) * 3.4562850092657564), |z| ≈ 5 > 3
As expected, we lose many digits precision due to massive cancellation.

...

>>> CiS(3.9, erf = lambda x: 1 - erfc(x))
(0.42233271026093344 + 0.47520240235068834j)

Bringing up an old thread with a question and a comment. The question: How do you get an accurate erfc(z) from 1-erf(z)? erf(3.9) = 0.999999965208, and that number subtracted from 1 is 0.000000034792 which has only 5 significant digits.

Now the comment: Much of the loss of precision in the first example comes from the large negative x in hyp1f1(1/2, 3/2, -x*x) which involves summing large alternating positive and negative numbers resulting in cancellation losses. The result, 0.22723763473 has only 7 correct digits. J-M Baillard's Special Function package uses a different formula, erf x = (2x/Pi^(1/2)) * exp(-x*x) * 1F1( 1 , 3/2 , x^2 ), which returns 0.999999965204, which is within 4 ULPs of the exact result. The positive x seems to result in more stable behavior of the hypergeometric 1F1 program.

Additional note: Baillard's program for C(x) and S(x) also uses a continued fraction if x > sqrt(pi).

[Albert Chan]

(04-04-2024 08:35 PM)John Keith Wrote:  How do you get an accurate erfc(z) from 1-erf(z)? erf(3.9) = 0.999999965208,
and that number subtracted from 1 is 0.000000034792 which has only 5 significant digits.

For big z, it is the other way around, erf(z) = 1 - erfc(z)
We need only few sig. digits of erfc(z) to get accurate erf(z)

Quote:erf x = (2x/Pi^(1/2)) * exp(-x*x) * 1F1( 1 , 3/2 , x^2 )

For CiS(3.9), erf argument is complex, A&S eqn 7.1.6 may not help with accuracy.

p2> k = 2/sqrt(pi)
p2> erf = lambda x: k*x * exp(-x*x) * hyp1f1(1, 3/2, x*x)
p2> CiS = lambda x, erf=erf: (1+1j)/2 * erf((1-1j)/k * x)
p2> CiS(3.9)
(0.42233255931844371+0.47520226114155328j)

[John Keith]

(04-04-2024 10:53 PM)Albert Chan Wrote:  For CiS(3.9), erf argument is complex, A&S eqn 7.1.6 may not help with accuracy.

Thanks, that makes sense. I guess that's why continued fraction evaluation is needed for large z.

[Albert Chan]

FYI, there is a problem with erfc continued fraction formula.
Calculate bottom up, if any CF denominators caused massive cancellations, all we get is garbage. (*)

Cas> erf(z := 3.16845945394*i)

1 - 58608047158.3*i     // erf(z) = 1 - erfc(z) ... but wrong

Cas> 2/sqrt(pi) * quad(t -> exp(-t*t), 0, z)[1]

4325.2249097*i            // actual integration, erf(z) correct.

(11-16-2023 01:22 PM)parisse Wrote:  I have disabled the continued fraction, it should now work correctly
(it will just be a little bit slower if |z| is between 4 and 6.5).

(*) With re(z)≥0, denominator may get to 0 only if re(z)=0

z + n/z = r*cis(θ) + n/r*cis(-θ)

Since cos(-θ) = cos(θ), we have:
--> re(z + n/z) ≥ re(z)
--> re(z + (n-1/2)/(z + n/z)) ≥ re(z)
...
--> re(z + (1/2)/(z + 1/(z + (3/2)/(z + ...)))) ≥ re(z)