Riemann's Zeta Function - another approach (RPL)

[Dieter]

(07-25-2017 04:00 PM)Gerson W. Barbosa Wrote:  If it is just a copy & paste matter, would you please provide a listing? Thanks!

I am currently experimenting with some adjustments to squeeze out the best possible accuracy. Here is the latest version. It includes several measures to keep intermediate results with more digits than required by the final result. That's why the constant c0 in line 35 and 119 has been decreased by 0,57 (which is added back later) so that in effect this crucial value can be given to 12 decimals.

Here is the listing:

Code:

 01 LBL "ZETA"
 02 STO 00
 03 SQRT
 04 RCL 00
 05 1
 06 -
 07 1/X
 08 LASTX
 09 X<0?
 10 GTO 97
 11 2
 12 RCL 00
 13 X>Y?
 14 GTO 96
 15 LASTX
 16 LASTX
 17 LASTX
 18 -1.276 E-8
 19 *
 20 7.05133 E-6
 21 -
 22 *
 23 9.721157 E-5
 24 +
 25 *
 26 3.4243368 E-4
 27 -
 28 *
 29 4.84515482 E-3
 30 -
 31 *
 32 7.281584288 E-2
 33 +
 34 *
 35 7.215664988 E-3
 36 +
 37 GTO 98
 38 LBL 96
 39 24
 40 RCL 00
 41 /
 42 2
 43 +
 44 INT
 45 ST+ X
 46 22
 47 X>Y?
 48 X<>Y
 49 STO 01
 50 RCL 00
 51 CHS
 52 STO 00
 53 CLX
 54 LBL 01
 55 RCL Y
 56 RCL 00
 57 Y^X
 58 -
 59 CHS
 60 DSE Y
 61 GTO 01
 62 RCL 00
 63 ST+ X
 64 1
 65 -
 66 RCL 01
 67 X^2
 68 24
 69 *
 70 /
 71 1
 72 RCL 00
 73 -
 74 8
 75 /
 76 RCL 01
 77 /
 78 +
 79 .5
 80 +
 81 RCL 01
 82 +
 83 RCL 00
 84 Y^X
 85 2
 86 /
 87 +
 88 1
 89 RCL 00
 90 +
 91 2
 92 LN
 93 *
 94 E^X-1
 95 CHS
 96 /
 97 GTO 99
 98 LBL 97
 99 ENTER
100 ENTER
101 ENTER
102 -8.4715 E-7
103 *
104 7.51334 E-6
105 -
106 *
107 9.609657 E-5
108 +
109 *
110 3.42683396 E-4
111 -
112 *
113 4.84527616 E-3
114 -
115 *
116 7.281583446 E-2
117 +
118 *
119 7.215664464 E-3
120 +
121 LBL 98
122 RDN
123 1/X
124 INT
125 ST* Z
126 SIGN
127 ST/ X
128 ST- Z
129 X<> L
130 RDN
131 /
132 -
133 R^
134 .57
135 +
136 +
137 LBL 99
138 END

The program uses a polynomial approximation for 0≤x<1, another one for 1<x≤2 and your original method for x>2. Line 03 and 07 generate an error if x<0 or x=1. The steps between LBL 98 and 99 add 1/(x–1) + 0,57 while trying to preseve as much accuracy as possible. This way e.g. Zeta(1,5) carries 12 digits (...7534869) before the last addition delivers the final result 2,612375349.

Dieter

[Gerson W. Barbosa]

(07-25-2017 04:29 PM)Dieter Wrote:  Here is the listing:

Code:

 01 LBL "ZETA"
 02 STO 00
 03 SQRT
 04 RCL 00
 05 1
 06 -
 07 1/X
 08 LASTX
 09 X<0?
 10 GTO 97
 11 2
 12 RCL 00
 13 X>Y?
 14 GTO 96
 15 LASTX
 16 LASTX
 17 LASTX
 18 -1.276 E-8
 19 *
 20 7.05133 E-6
 21 -
 22 *
 23 9.721157 E-5
 24 +
 25 *
 26 3.4243368 E-4
 27 -
 28 *
 29 4.84515482 E-3
 30 -
 31 *
 32 7.281584288 E-2
 33 +
 34 *
 35 7.215664988 E-3
 36 +
 37 GTO 98
 38 LBL 96
 39 24
 40 RCL 00
 41 /
 42 2
 43 +
 44 INT
 45 ST+ X
 46 22
 47 X>Y?
 48 X<>Y
 49 STO 01
 50 RCL 00
 51 CHS
 52 STO 00
 53 CLX
 54 LBL 01
 55 RCL Y
 56 RCL 00
 57 Y^X
 58 -
 59 CHS
 60 DSE Y
 61 GTO 01
 62 RCL 00
 63 ST+ X
 64 1
 65 -
 66 RCL 01
 67 X^2
 68 24
 69 *
 70 /
 71 1
 72 RCL 00
 73 -
 74 8
 75 /
 76 RCL 01
 77 /
 78 +
 79 .5
 80 +
 81 RCL 01
 82 +
 83 RCL 00
 84 Y^X
 85 2
 86 /
 87 +
 88 1
 89 RCL 00
 90 +
 91 2
 92 LN
 93 *
 94 E^X-1
 95 CHS
 96 /
 97 GTO 99
 98 LBL 97
 99 ENTER
100 ENTER
101 ENTER
102 -8.4715 E-7
103 *
104 7.51334 E-6
105 -
106 *
107 9.609657 E-5
108 +
109 *
110 3.42683396 E-4
111 -
112 *
113 4.84527616 E-3
114 -
115 *
116 7.281583446 E-2
117 +
118 *
119 7.215664464 E-3
120 +
121 LBL 98
122 RDN
123 1/X
124 INT
125 ST* Z
126 SIGN
127 ST/ X
128 ST- Z
129 X<> L
130 RDN
131 /
132 -
133 R^
134 .57
135 +
136 +
137 LBL 99
138 END

Very nice! Now everything can be done in less than 30 seconds!

Code:


01 LBL "Z"
02 X<0?
03 GTO 01
04 XEQ "ZETA"
05 GTO 02
06 LBL 01
07 1
08 -
09 PI
10 ST+ X
11 X<>Y
12 Y^X
13 STO 05
14 LASTX
15 CHS
16 XEQ "ZETA"
17 ST* 05
18 RCL 00
19 ABS
20 1
21 -
22 XEQ "GAM+1"
23 ST* 05
24 RCL 00
25 CHS
26 1
27 ASIN
28 *
29 SIN
30 RCL 05
31 *
32 ST+ X
33 LBL 02
34 END

"ZETA" in lines 04 and 16 is your program above. "GAM+1" in line 22 is your program here.

 3 XEQ "Z"  -->   1.202056903         (15.7 s)
 2.001 R/S  -->   1.643997513     (2) (17.4 s)
 2     R/S  -->   1.644934067         ( 4.7 s)
 1.5   R/S  -->   2.612375349         ( 4.8 s)
 0.5   R/S  -->  -1.460354509         ( 5.4 s)
 0     R/S  -->  -0.500000000         ( 4.7 s)
-0.5   R/S  -->  -0.2078862248   (50) (15.4 s)
-1     R/S  -->  -0.08333333322  (33) (14.6 s)
-1.001 R/S  -->  -0.08316803706  (23) (27.4 s)
-1.5   R/S  -->  -0.02548520184  (90) (27.3 s)
-2     R/S  -->   0.00000000000       (28.8 s)
-3     R/S  -->   0.008333333324 (33) (23.3 s)
-5     R/S  -->  -0.003968253964  (8) (21.2 s)
-7     R/S  -->   0.004166666670 (67) (20.1 s)

[Dieter]

(07-29-2017 03:37 PM)Gerson W. Barbosa Wrote:  "ZETA" in lines 04 and 16 is your program above. "GAM+1" in line 22 is your program here.

Ah yes, the Gamma thread. There were some more promising approaches with even better accuracy than GAM+1. ;-)

BTW, I see your program has a RCL 00 in line 24. For x≤2 the ZETA routine leaves x in R00, but for x>2 R00 finally holds –x. Have you considered this?
Anyway, since the results look fine for smaller and larger negative x the program must be working fine. ;-)

BTW2, on 10-digit calculators 2 pi does not round very well, the last digit is 1 unit high. So I suggest to replace step 09/10 in your program with 360 D–R which returns the correct value 6,283185307. The same is true for pi/4 or pi/6 where 45 resp. 30 D–R yields ten correct digits.

Dieter

[Namir]

Good work Dieter!! I will type in your code and play with it.

Namir

[Gerson W. Barbosa]

(07-30-2017 07:57 AM)Dieter Wrote:  BTW, I see your program has a RCL 00 in line 24. For x≤2 the ZETA routine leaves x in R00, but for x>2 R00 finally holds –x. Have you considered this?

Yes, that's what ABS in line 19 is for.

(07-30-2017 07:57 AM)Dieter Wrote:  BTW2, on 10-digit calculators 2 pi does not round very well, the last digit is 1 unit high. So I suggest to replace step 09/10 in your program with 360 D–R which returns the correct value 6,283185307. The same is true for pi/4 or pi/6 where 45 resp. 30 D–R yields ten correct digits.

That's a good suggestion, but we'd need x! (or Gamma) to be that accurate too. Is there a math module with x! or Gamma?

A few guard digits (perhaps just a couple of them) combined with built-in Gamma might give perfect 10-digit results most always, even when using 10-digits constants, which is quite impressive. On Free42:

 3 XEQ "Z"  -->   1.20205690313 (6)         
 2.001 R/S  -->   1.64399751259 (24)
 2     R/S  -->   1.64493406685         
 1.5   R/S  -->   2.61237534868 (9)
 0.5   R/S  -->  -1.46035450879 (81)
 0     R/S  -->  -0.50000000000
-0.5   R/S  -->  -0.207886224977
-1     R/S  -->  -0.0833333333333
-1.001 R/S  -->  -0.0831680372461 (281)
-1.5   R/S  -->  -0.0254852018937 (898)
-2     R/S  -->   0.0000000000000       
-3     R/S  -->   0.00833333333191 (333)
-5     R/S  -->  -0.00396825396801 (25)
-7     R/S  -->   0.00416666666663 (7)

Code:


00 { 57-Byte Prgm }
01▸LBL "Z"
02 X<0?
03 GTO 01
04 XEQ "ZETA"
05 GTO 02
06▸LBL 01
07 1
08 -
09 PI
10 STO+ ST X
11 X<>Y
12 Y↑X
13 STO 05
14 LASTX
15 +/-
16 GAMMA
17 STO× 05
18 LASTX
19 XEQ "ZETA"
20 STO× 05
21 1
22 RCL 00
23 ABS
24 -
25 1
26 ASIN
27 ×
28 SIN
29 RCL× 05
30 STO+ ST X
31▸LBL 02
32 END

Code:



00 { 359-Byte Prgm }
01▸LBL "ZETA"
02 STO 00
03 SQRT
04 -1
05 RCL+ 00
06 1/X
07 LASTX
08 X<0?
09 GTO 97
10 2
11 RCL 00
12 X>Y?
13 GTO 96
14 LASTX
15 LASTX
16 LASTX
17 -1,276ᴇ-8
18 ×
19 7,05133ᴇ-6
20 -
21 ×
22 9,721157ᴇ-5
23 +
24 ×
25 3,4243368ᴇ-4
26 -
27 ×
28 0,00484515482
29 -
30 ×
31 0,07281584288
32 +
33 ×
34 0,007215664988
35 +
36 GTO 98
37▸LBL 96
38 24
39 RCL 00
40 ÷
41 2
42 +
43 IP
44 STO+ ST X
45 22
46 X>Y?
47 X<>Y
48 STO 01
49 RCL 00
50 +/-
51 STO 00
52 CLX
53▸LBL 01
54 RCL ST Y
55 RCL 00
56 Y↑X
57 -
58 +/-
59 DSE ST Y
60 GTO 01
61 RCL 00
62 STO+ ST X
63 1
64 -
65 RCL 01
66 X↑2
67 24
68 ×
69 ÷
70 1
71 RCL- 00
72 8
73 ÷
74 RCL÷ 01
75 +
76 0,5
77 +
78 RCL+ 01
79 RCL 00
80 Y↑X
81 2
82 ÷
83 +
84 1
85 RCL+ 00
86 2
87 LN
88 ×
89 E↑X-1
90 +/-
91 ÷
92 GTO 99
93▸LBL 97
94 ENTER
95 ENTER
96 ENTER
97 -8,4715ᴇ-7
98 ×
99 7,51334ᴇ-6
100 -
101 ×
102 9,609657ᴇ-5
103 +
104 ×
105 3,42683396ᴇ-4
106 -
107 ×
108 0,00484527616
109 -
110 ×
111 0,07281583446
112 +
113 ×
114 0,007215664464
115 +
116▸LBL 98
117 R↓
118 1/X
119 IP
120 STO× ST Z
121 SIGN
122 STO÷ ST X
123 STO- ST Z
124 X<> ST L
125 R↓
126 ÷
127 -
128 R↑
129 0,57
130 +
131 +
132▸LBL 99
133 END

[Gerson W. Barbosa]

(07-30-2017 04:17 PM)Gerson W. Barbosa Wrote:  

Code:


00 { 57-Byte Prgm }
01▸LBL "Z"
02 X<0?
03 GTO 01
04 XEQ "ZETA"
05 GTO 02
06▸LBL 01
07 1
08 -
09 PI
10 STO+ ST X
11 X<>Y
12 Y↑X
13 STO 05
14 LASTX
15 +/-
16 GAMMA
17 STO× 05
18 LASTX
19 XEQ "ZETA"
20 STO× 05
21 1
22 RCL 00
23 ABS
24 -
25 1
26 ASIN
27 ×
28 SIN
29 RCL× 05
30 STO+ ST X
31▸LBL 02
32 END

Or, slightly less inelegant,

Code:


00 { 55-Byte Prgm }
01▸LBL "Z"
02 X<0?
03 GTO 01
04 XEQ "ZETA"
05 GTO 02
06▸LBL 01
07 2
08 X<>Y
09 Y↑X
10 PI
11 LASTX
12 1
13 -
14 Y↑X
15 LASTX
16 +/-
17 GAMMA
18 LASTX
19 R↓
20 ×
21 ×
22 1
23 ASIN
24 RCL× ST Z
25 COS
26 ×
27 STO 05
28 R↑
29 XEQ "ZETA"
30 RCL× 05
31▸LBL 02
32 END

[/quote]

[Dieter]

(07-30-2017 04:17 PM)Gerson W. Barbosa Wrote:  

(07-30-2017 07:57 AM)Dieter Wrote:  BTW, I see your program has a RCL 00 in line 24. For x≤2 the ZETA routine leaves x in R00, but for x>2 R00 finally holds –x. Have you considered this?

Yes, that's what ABS in line 19 is for.

I see, but the HP-41 version has is another RCL 00 without an ABS in line 24. ;-)

(07-30-2017 04:17 PM)Gerson W. Barbosa Wrote:  That's a good suggestion, but we'd need x! (or Gamma) to be that accurate too. Is there a math module with x! or Gamma?

I do not know of an (official HP)-ROM with a full-precision 10-digit Mcode-Gamma-implementation.

(07-30-2017 04:17 PM)Gerson W. Barbosa Wrote:  A few guard digits (perhaps just a couple of them) combined with built-in Gamma might give perfect 10-digit results most always, even when using 10-digits constants, which is quite impressive.

To assess the final accuracy it might be helpful to know the error of the two polynomial approxmations for x between 0 and 2. If evaluated with sufficient precision and using the coefficients given in the HP-41 program (note that c0 effectively has 12 digits since 0,57 is added later) the one for 0≤x<1 has a largest error of ~3,5 units in the 12th significant digit, while for 1<x≤2 it's less than 0,7 units.

(07-30-2017 04:17 PM)Gerson W. Barbosa Wrote:  On Free42:

I tried Free42 BCD where Gamma seems to be good for 30+ digits. So the results should only be limited by the approximation error. However, for x>2 the number of terms (cf. line 38 ff. in the 42s Zeta program) of course has to be adjusted in a higher precision environment. ;-) This should improve the results for x>2 resp. x<–1.

Dieter

[Gerson W. Barbosa]

(07-30-2017 10:18 PM)Dieter Wrote:  

(07-30-2017 04:17 PM)Gerson W. Barbosa Wrote:  Yes, that's what ABS in line 19 is for.

I see, but the HP-41 version has is another RCL 00 without an ABS in line 24. ;-)

Yes, but it's located after a call to your "GAM+1" program, which incidentally saves the argument in register 0 also. But you're not supposed to have to remember that after 3+ years. I'd failed to notice what your concern was about, sorry!

(07-30-2017 10:18 PM)Dieter Wrote:  

(07-30-2017 04:17 PM)Gerson W. Barbosa Wrote:  A few guard digits (perhaps just a couple of them) combined with built-in Gamma might give perfect 10-digit results most always, even when using 10-digits constants, which is quite impressive.

To assess the final accuracy it might be helpful to know the error of the two polynomial approxmations for x between 0 and 2. If evaluated with sufficient precision and using the coefficients given in the HP-41 program (note that c0 effectively has 12 digits since 0,57 is added later) the one for 0≤x<1 has a largest error of ~3,5 units in the 12th significant digit, while for 1<x≤2 it's less than 0,7 units.

(07-30-2017 04:17 PM)Gerson W. Barbosa Wrote:  On Free42:

I tried Free42 BCD where Gamma seems to be good for 30+ digits. So the results should only be limited by the approximation error. However, for x>2 the number of terms (cf. line 38 ff. in the 42s Zeta program) of course has to be adjusted in a higher precision environment. ;-) This should improve the results for x>2 resp. x<–1.

This was just a quick test on Free42. I have yet to try it on the 42S, where some accuracy loss is expected, but not so much, I hope. Yes, the number of terms will have to be increased accordingly. Despite that, the running times will be better, since the 42S is about five times faster than the 41.

Gerson.

[Gerson W. Barbosa]

(07-30-2017 10:18 PM)Dieter Wrote:  I tried Free42 BCD where Gamma seems to be good for 30+ digits. So the results should only be limited by the approximation error. However, for x>2 the number of terms (cf. line 38 ff. in the 42s Zeta program) of course has to be adjusted in a higher precision environment. ;-) This should improve the results for x>2 resp. x<–1.

I've tested it on Emu42 set to "Authentic Calculator Speed". The number of terms hasn't been adjusted yet, thus accuracy differences between both can be highlighted:


 3 XEQ "ZETA" -->   1.20205690313      (6) (4.0 s)
 2.001 R/S    -->   1.64399751259     (24) (4.4 s)
 2     R/S    -->   1.64493406684      (5) (1.2 s)
 1.5   R/S    -->   2.61237534867      (9) (1.2 s)
 0.5   R/S    -->  -1.46035450879     (81) (1.1 s)
 0     R/S    -->  -0.500000000004     (0) (1.1 s)
-0.5   R/S    -->  -0.207886224976     (7) (1.6 s)
-1     R/S    -->  -0.0833333333333   (33) (1.5 s)
-1.001 R/S    -->  -0.0831680372464  (281) (4.8 s)
-1.5   R/S    -->  -0.0254852018936  (898) (4.7 s)
-2     R/S    -->   0.0000000000000        (4.4 s)
-3     R/S    -->   0.00833333333194 (333) (3.7 s)
-5     R/S    -->  -0.00396825396798 (825) (3.3 s)
-7     R/S    -->   0.00416666666664   (7) (3.1 s)

Code:


00 { 402-Byte Prgm }
01▸LBL "ZETA"
02 X≥0?
03 GTO 00
04 1
05 -
06 PI
07 STO+ ST X
08 X<>Y
09 Y↑X
10 LASTX
11 +/-
12 GAMMA
13 LASTX
14 R↓
15 ×
16 1
17 ASIN
18 RCL× ST T
19 COS
20 ×
21 STO 02
22 R↑
23 XEQ 00
24 RCL× 02
25 STO+ ST X
26 GTO 99
27▸LBL 00
28 STO 00
29 SQRT
30 -1
31 RCL+ 00
32 1/X
33 LASTX
34 X<0?
35 GTO 97
36 2
37 RCL 00
38 X>Y?
39 GTO 96
40 LASTX
41 LASTX
42 LASTX
43 -1.276ᴇ-8
44 ×
45 7.05133ᴇ-6
46 -
47 ×
48 9.721157ᴇ-5
49 +
50 ×
51 3.4243368ᴇ-4
52 -
53 ×
54 0.00484515482
55 -
56 ×
57 0.07281584288
58 +
59 ×
60 0.007215664988
61 +
62 GTO 98
63▸LBL 96
64 24
65 RCL 00
66 ÷
67 2
68 +
69 IP
70 STO+ ST X
71 22
72 X>Y?
73 X<>Y
74 STO 01
75 RCL 00
76 +/-
77 STO 00
78 CLX
79▸LBL 95
80 RCL ST Y
81 RCL 00
82 Y↑X
83 -
84 +/-
85 DSE ST Y
86 GTO 95
87 RCL 00
88 STO+ ST X
89 1
90 -
91 RCL 01
92 X↑2
93 24
94 ×
95 ÷
96 1
97 RCL- 00
98 8
99 ÷
100 RCL÷ 01
101 +
102 0.5
103 +
104 RCL+ 01
105 RCL 00
106 Y↑X
107 2
108 ÷
109 +
110 1
111 RCL+ 00
112 2
113 LN
114 ×
115 E↑X-1
116 +/-
117 ÷
118 GTO 99
119▸LBL 97
120 ENTER
121 ENTER
122 ENTER
123 -8.4715ᴇ-7
124 ×
125 7.51334ᴇ-6
126 -
127 ×
128 9.609657ᴇ-5
129 +
130 ×
131 3.42683396ᴇ-4
132 -
133 ×
134 0.00484527616
135 -
136 ×
137 0.07281583446
138 +
139 ×
140 0.007215664464
141 +
142▸LBL 98
143 R↓
144 1/X
145 IP
146 STO× ST Z
147 SIGN
148 STO÷ ST X
149 STO- ST Z
150 X<> ST L
151 R↓
152 ÷
153 -
154 R↑
155 0.57
156 +
157 +
158▸LBL 99
159 END

Your code (lines 027-159) is still intact :-)

Gerson.

[Dieter]

(07-31-2017 03:59 PM)Gerson W. Barbosa Wrote:  I've tested it on Emu42 set to "Authentic Calculator Speed". The number of terms hasn't been adjusted yet, thus accuracy differences between both can be highlighted:

Hmmm... x=0 should return exactly –0,5 because the coefficients are designed accordingly. So I wonder why you get –0,500000000004. This would mean that one of the coefficients is off by 4 E-12. ?!?

BTW the optimized coefficient set has two slight changes in the last digit:
line 123: -8.47149 E-7
line 131: 3.42683395 E-4

And the error check for x<0 (SQRT in line 29) can be omitted here. ;-)

But remember, the two polynomial approximations are intended for 10 digit calculators. While the approximation for 1<x≤2 has an error less than 1 unit in the 12th digit and thus is useable for the 42s, the one for 0≤x<1 may be off by 4 ULP, so here one more term would make sense. Looking at your examples, indeed the largest error between –1 and +2 is 2 ULP.

Dieter

[Gerson W. Barbosa]

(07-31-2017 05:00 PM)Dieter Wrote:  

(07-31-2017 03:59 PM)Gerson W. Barbosa Wrote:  I've tested it on Emu42 set to "Authentic Calculator Speed". The number of terms hasn't been adjusted yet, thus accuracy differences between both can be highlighted:

Hmmm... x=0 should return exactly –0,5 because the coefficients are designed accordingly. So I wonder why you get –0,500000000004. This would mean that one of the coefficients is off by 4 E-12. ?!?

BTW the optimized coefficient set has two slight changes in the last digit:
line 123: -8.47149 E-7
line 131: 3.42683395 E-4

I wonder what might have caused these differences. The coefficients were copied from your HP-41C listing:

 43 -1.276ᴇ-8            18 -1.276 E-8
 44 ×                    19 *
 45 7.05133ᴇ-6           20 7.05133 E-6
 46 -                    21 -
 47 ×                    22 *
 48 9.721157ᴇ-5          23 9.721157 E-5
 49 +                    24 +
 50 ×                    25 *
 51 3.4243368ᴇ-4         26 3.4243368 E-4
 52 -                    27 -
 53 ×                    28 *
 54 0.00484515482        29 4.84515482 E-3
 55 -                    30 -
 56 ×                    31 *
 57 0.07281584288        32 7.281584288 E-2
 58 +                    33 +
 59 ×                    34 *
 60 0.007215664988       35 7.215664988 E-3
 61 +                    36 +


123 -8.4715ᴇ-7          102 -8.4715 E-7
124 ×                   103 *
125 7.51334ᴇ-6          104 7.51334 E-6
126 -                   105 -
127 ×                   106 *
128 9.609657ᴇ-5         107 9.609657 E-5
129 +                   108 +
130 ×                   109 *
131 3.42683396ᴇ-4       110 3.42683396 E-4
132 -                   111 -
133 ×                   112 *
134 0.00484527616       113 4.84527616 E-3
135 -                   114 -
136 ×                   115 *
137 0.07281583446       116 7.281583446 E-2
138 +                   117 +
139 ×                   118 *
140 0.007215664464      119 7.215664464 E-3
141 +                   120 +

(07-31-2017 05:00 PM)Dieter Wrote:  And the error check for x<0 (SQRT in line 29) can be omitted here. ;-)

That was a blind copy, as you can see.

(07-31-2017 05:00 PM)Dieter Wrote:  But remember, the two polynomial approximations are intended for 10 digit calculators. While the approximation for 1<x≤2 has an error less than 1 unit in the 12th digit and thus is useable for the 42s, the one for 0≤x<1 may be off by 4 ULP, so here one more term would make sense. Looking at your examples, indeed the largest error between –1 and +2 is 2 ULP.

This is good enough for me, at least for the time being. I will only double the number of terms, which is an easy thing to do, and see what I get. I'll leave the 42S as it is so anyone interested can improve it, starting by deleting line 29 :-)

Thanks for your work on the two polynomial approximations which have made this fast enough for all arguments in the valid range.

Gerson.

[Dieter]

(07-31-2017 06:31 PM)Gerson W. Barbosa Wrote:  I wonder what might have caused these differences. The coefficients were copied from your HP-41C listing:

Could you check what the polynomial returns for x=0 at the point directly before LBL 98 ? Here the value should be exactly –0,07.

But... since you use the reflection formula the approximation only has to cover the range from x=0,5 to 1. This can be done with one power less:

Zeta(x) ~ 1/u + 0,577215664944 + 0,07281584734 u – 0,00484515533 u² – 0,00034213951 u³ + 0,00009740462 u⁴ – 0,000005872685 u⁵
where u = x – 1  and  0,5 ≤ x < 1.

Within the given domain this approximation has a largest error of ~0,3 units in the 12th significant digit.

If you replace the 0...1 polynomial with this one be sure to set the last coefficient (constant term) to 0,007215664944 since the missing 0,57 is added later.
On a 12-digit calculator you can even add two more digits at make it ...494433 while the u⁵ coefficient becomes ...72675. This way the error practically matches the theoretical optimum at ~0,27 ULP.

Dieter

[Dieter]

(07-31-2017 06:31 PM)Gerson W. Barbosa Wrote:  I wonder what might have caused these differences. The coefficients were copied from your HP-41C listing:

...
 60 0.007215664988       35 7.215664988 E-3
...
140 0.007215664464      119 7.215664464 E-3

I just found that this might be the reason why your result for x=0 is not exactly –0,5 and others may be slightly off as well: for instance Zeta(0,5) and Zeta(2) should return the exact 12-digit values -1.46035450881 and 1,64493406685, actually Free42 returns even 13 correct digits here.

On Free42 (and I suppose on a real hardware 42s) you cannot enter line 60 resp. 140 as the numbers are too long and the last digit is omitted. Could you check this please and try 7.215664988 E-3 resp. 7.215664464 E-3 here?

BTW, I am now testing a Free42 version, essentially your code but with a new approximation for 0,5...1. The only problem is x=0 as the reflection formula would cause an attempt at calculating Zeta(1). Anyway, the results look quite promising once the number of terms for x>2 is sufficiently large. With about 200/x terms (a very rough estimate) the largest error I noticed for your examples above was 2 ULP.

Dieter

[Gerson W. Barbosa]

(08-01-2017 11:46 AM)Dieter Wrote:  

(07-31-2017 06:31 PM)Gerson W. Barbosa Wrote:  I wonder what might have caused these differences. The coefficients were copied from your HP-41C listing:

...
 60 0.007215664988       35 7.215664988 E-3
...
140 0.007215664464      119 7.215664464 E-3

I just found that this might be the reason why your result for x=0 is not exactly –0,5 and others may be slightly off as well: for instance Zeta(0,5) and Zeta(2) should return the exact 12-digit values -1.46035450881 and 1,64493406685, actually Free42 returns even 13 correct digits here.

On Free42 (and I suppose on a real hardware 42s) you cannot enter line 60 resp. 140 as the numbers are too long and the last digit is omitted. Could you check this please and try 7.215664988 E-3 resp. 7.215664464 E-3 here?

Bingo! Now I get -0.5 exactly. Actually, Free 42 accepts numbers with more than 30 digits.

Have to leave now. Will comment later.

Gerson.

PS: Here is the table after the suggested modifications:

 2     R/S    -->   1.64493406685      
 1.5   R/S    -->   2.61237534868      (9) 
 0.5   R/S    -->  -1.46035450879     (81)
 0     R/S    -->  -0.500000000000
-0.5   R/S    -->  -0.207886224978     (7) 
-1     R/S    -->  -0.0833333333334    (3)

  60  7.215664988E-3 
 123 -8.47149E-7
 131  3.42683395E-4
 140  7.215664464E-3

That was difficult to detect because Emu48 shows the lines


 60 0.007215664988
140 0.007215664464

Yet they could not have been entered on the real 42S.
The .RAW file I used was exported by Free42.

From now on I'll take care when exporting Free42 .RAW files to Emu42.

Code:


00 { 45-Byte Prgm }
01▸LBL "T"
02 0.333333333333333333333333333333333
03 1/X
04 3
05 -
06 END

XEQ "T" --> 3.E-33

[Gerson W. Barbosa]

(08-01-2017 11:46 AM)Dieter Wrote:  BTW, I am now testing a Free42 version, essentially your code but with a new approximation for 0,5...1. The only problem is x=0 as the reflection formula would cause an attempt at calculating Zeta(1).

This is the formula I used in the previous program, but then only when x<0:

\[\zeta(x)=2\cdot\left ( 2\pi \right )^{-(1-x)}\cdot \cos\left ( (1-x) \sin^{-1}(1)\right )\cdot \Gamma (1-x)\cdot \zeta (1-x)\]


I remember I had to hard-code Zeta(0) = -1/2 here (line 116).

Gerson.

[Dieter]

(08-01-2017 04:54 PM)Gerson W. Barbosa Wrote:  I remember I had to hard-code Zeta(0) = -1/2 here (line 116).

The problem is any x sufficiently close to zero as 1–x may round to 1. This also seems to be a problem on the WP34s, cf. my post in the other forum. Here the polynomial approximation that goes down to x=0 has an advantage.

In the meantime I have looked at the required number of terms for the iterative method (x>2). With about 400/x² + 10 terms the remaining error is roughly 2 units in the 13th digit which is as good as the best of the polynomial approximations. But this requires more than 100 terms. Then I found another quite effective method that requires less terms: n = 2*int(100/x^1,5 + 5) uses merely 80 terms at x=2 and has an error that can be easily compensated with a simple formula: for x<4 subtract 4,4 E–12*(4–x) from the final result. This should limit the error to less than ±1 unit in the 13th digit. And the correction can be done in just 7 steps.

This can also be done for 10-digit accuracy, e.g. the HP-41 version. Here one may use n = 2*int(22/x + 4) terms and decrease the result for x<3 by 5 E–10*(3–x). This should limit the error to about one unit in the 11th digit, mostly even less. At least that's what can be achieved with sufficient arithmetic precision. On a physical HP-41 or a true emulator this does not make any difference as the correction is always less than 5 E–10 in a result >1. ;-) But I will change the number of terms to 2*int(24/x+3) and remove the limitation to n=22.

Dieter

[Gerson W. Barbosa]

(08-02-2017 11:19 AM)Dieter Wrote:  

(08-01-2017 04:54 PM)Gerson W. Barbosa Wrote:  I remember I had to hard-code Zeta(0) = -1/2 here (line 116).

The problem is any x sufficiently close to zero as 1–x may round to 1. This also seems to be a problem on the WP34s, cf. my post in the other forum. Here the polynomial approximation that goes down to x=0 has an advantage.

It appears the problem has been fixed on the WP34S. Without your testings, probably it would remain unnoticed for a long time.

Back to the HP-41:

Code:



 01▸LBL "ZETA"
 02 X>0?
 03 GTO 00
 04 X=0?
 05 GTO 00
 06 CHS
 07 1
 08 +
 09 STO 02
 10 XEQ "GAMMA"
 11 STO 03
 12 RCL 02
 13 XEQ 00
 14 RCL 03
 15 *
 16 PI
 17 ST+ X
 18 RCL 02
 19 Y^X
 20 /
 21 1
 22 ASIN
 23 RCL 02
 24 *
 25 COS
 26 *
 27 ST+ X
 28 GTO 99
 29▸LBL 00
 30 STO 00
 31 -1
 32 +
 33 1/X
 34 LASTX
 35 X<0?
 36 GTO 97
 37 2
 38 RCL 00
 39 X>Y?
 40 GTO 96
 41 LASTX
 42 LASTX
 43 LASTX
 44 -1.276E -8
 45 *
 46 7.05133E -6
 47 -
 48 *
 49 9.721157E -5
 50 +
 51 *
 52 3.4243368E -4
 53 -
 54 *
 55 4.84515482E -3
 56 -
 57 *
 58 7.281584288E -2
 59 +
 60 *
 61 7.21566498E -3
 62 +
 63 GTO 98
 64▸LBL 96
 65 24
 66 RCL 00
 67 /
 68 2
 69 +
 70 INT
 71 ST+ X
 72 22
 73 X>Y?
 74 X<>Y
 75 STO 01
 76 RCL 00
 77 CHS
 78 STO 00
 79 CLX
 80▸LBL 01
 81 RCL Y
 82 RCL 00
 83 Y^X
 84 -
 85 CHS
 86 DSE Y
 87 GTO 01
 88 RCL 00
 89 ST+ X
 90 1
 91 -
 92 RCL 01
 93 X^2
 94 24
 95 *
 96 /
 97 1
 98 RCL 00
 99 -
100 8
101 /
102 RCL 01
103 /
104 +
105 .5
106 +
107 RCL 01
108 +
109 RCL 00
110 Y^X
111 2
112 /
113 +
114 1
115 RCL 00
116 +
117 2
118 LN
119 *
120 E^X-1
121 CHS
122 /
123 GTO 99
124▸LBL 97
125 ENTER^
126 ENTER^
127 ENTER^
128 -8.4715E -7
129 *
130 7.51334E -6
131 -
132 *
133 9.609657E -5
134 +
135 *
136 3.42683396E -4
137 -
138 *
139 4.84527616E -3
140 -
141 *
142 7.281583446E -2
143 +
144 *
145 7.215664464E -3
146 +
147▸LBL 98
148 RDN
149 1/X
150 INT
151 ST* Z
152 SIGN
153 STO/ X
154 STO- Z
155 X<> L
156 RDN
157 /
158 -
159 R^
160 .57
161 +
162 +
163▸LBL 99
164 END

Code:


 01▸LBL "GAMMA"
 02 -1
 03 X<>Y
 04 +
 05 1
 06 STO 00
 07 X>Y?
 08 GTO 01
 09 ST- L
 10 LASTX
 11 STO 00
 12 INT
 13 1
 14▸LBL 00
 15 RCL 00
 16 *
 17 1
 18 ST- 00
 19 RDN
 20 DSE Y
 21 GTO 00
 22 STO 00
 23▸LBL 01
 24 LASTX
 25 ENTER^
 26 ENTER^
 27 ENTER^
 28 1.604589926 E-2
 29 *
 30 2.641400819 E-1
 31 -
 32 *
 33 1.96580515
 34 +
 35 *
 36 8.729327662
 37 -
 38 *
 39 25.69590147
 40 +
 41 *
 42 52.63472435
 43 -
 44 *
 45 76.53826433
 46 +
 47 *
 48 78.92639573
 49 -
 50 *
 51 56.50761084
 52 +
 53 *
 54 26.37243835
 55 -
 56 *
 57 7.203398486
 58 +
 59 RCL 00
 60 *
 61 END

ZETA, from line 29 on, is essentially your code. GAMMA uses a 10th degree polynomial approximation for the interval [1..2] and simple multiplications for x > 2, but a higher order polynomial might be necessary for 10-digit accuracy. Anyway, it can be replaced with better GAMMA implementations. Still, the ideal 30-second limit is surpassed only for large negative arguments, close to the end of the valid range:

  3 XEQ ZETA  -->   1.202056903           (15.6 s)
  2.001 R/S   -->   1.643997513       (2) (17.0 s)
  2     R/S   -->   1.644934067           ( 4.5 s)
  1.5   R/S   -->   2.612375349           ( 4.5 s)
  0.5   R/S   -->  -1.460354509           ( 5.1 s)
  0     R/S   -->  -0.500000000           ( 4.6 s)
 -0.5   R/S   -->  -0.2078862450    (250) (12.8 s)
 -1     R/S   -->  -0.08333333384    (33) (11.6 s)
 -1.001 R/S   -->  -0.08316803696   (723) (24.3 s)
 -1.5   R/S   -->  -0.02548520436   (190) (24.0 s)
 -2     R/S   -->   0.00000000000         (23.6 s)
 -3     R/S   -->   0.008333333384   (33) (20.7 s)
 -5     R/S   -->  -0.003968253990   (68) (19.3 s)
 -7     R/S   -->   0.004166666686   (67) (18.4 s)
-15.16  R/S   -->   0.4964873534     (85) (18.5 s)
-33.34  R/S   -->  -1.924684098E10  (152) (21.9 s)
-41.42  R/S   -->  -3.506595602E16  (584) (24.1 s)
-48.49  R/S   -->  -3.653091058E22   (22) (26.2 s)
-58.59  R/S   -->   1.136304829E32  (792) (29.0 s)
-67.97  R/S   -->   1.832460467E40 (1182) (31.4 s)

[Dieter]

[quote='Gerson W. Barbosa' pid='76831' dateline='1501861065']
Back to the HP-41:

Code:

 01▸LBL "ZETA"
 02 X>0?
 03 GTO 00
 04 X=0?
 05 GTO 00
 ...

Just one remark: Testing condition1 OR condition2 can be done with a simple trick from the olden days. Simply have the inverse of condition1 followed by condition2.

Code:

 01▸LBL "ZETA"
 02 X≠0?
 03 X>0?
 04 GTO 00
 ...

And the calculation of the number of terms should be replaced by this:

Code:

 64▸LBL 96
 65 24
 66 RCL 00
 67 /
 68 4
 69 +
 70 INT
 71 ST+ X
 72 STO 01

At the moment I am trying a Free42 version. Due to the 34-digit precision some things can be substantially simplified, e.g. most of the code following LBL 98 can be replaced by a simple x<>y 1/x + without all the stack acrobatics. And a result which usually returns 12 valid digits should be possible.

Edit: here is the current version. It addresses the x~0 problem in a ...err... "creative" way in that the smallest |x| is set to 1E–20. #-)

Please note that this version is intended for high-precision environments like Free42 and may not perform as well on a regular 42s.

Code:

 00 { 416-Byte Prgm }
 01>LBL "ZETA"
 02 0.5
 03 X<>Y
 04 X>=Y?
 05 GTO 77
 06 ABS
 07 1E-20
 08 X<Y?
 09 X<>Y
 10 LASTX
 11 SIGN
 12 ×
 13 1
 14 X<>Y
 15 -
 16 XEQ 77
 17 PI
 18 STO+ ST X
 19 RCL 00
 20 Y^X
 21 ÷
 22 RCL 00
 23 GAMMA
 24 ×
 25 1
 26 ASIN
 27 RCL× 00
 28 COS
 29 ×
 30 STO+ ST X
 31 1
 32 RCL- 00
 33 STO 00
 34 Rv
 35 GTO 99
 36>LBL 77
 37 STO 00
 38 -1
 39 RCL+ 00
 40 1/X
 41 LASTX
 42 X<0?
 43 GTO 97
 44 2
 45 RCL 00
 46 X>Y?
 47 GTO 96
 48 LASTX
 49 LASTX
 50 LASTX
 51 -1.276E-8
 52 ×
 53 7.05133E-6
 54 -
 55 ×
 56 9.721157E-5
 57 +
 58 ×
 59 3.4243368E-4
 60 -
 61 ×
 62 0.00484515482
 63 -
 64 ×
 65 0.07281584288
 66 +
 67 ×
 68 7.215664988E-3
 69 +
 70 GTO 98
 71>LBL 96
 72 100
 73 RCL 00
 74 SQRT
 75 RCL× 00
 76 ÷
 77 5
 78 +
 79 IP
 80 STO+ ST X
 81 STO 01
 82 RCL 00
 83 +/-
 84 STO 00
 85 CLX
 86>LBL 95
 87 RCL ST Y
 88 RCL 00
 89 Y^X
 90 -
 91 +/-
 92 DSE ST Y
 93 GTO 95
 94 RCL 00
 95 STO+ ST X
 96 1
 97 -
 98 RCL 01
 99 X^2
100 24
101 ×
102 ÷
103 1
104 RCL- 00
105 8
106 ÷
107 RCL÷ 01
108 +
109 0.5
110 +
111 RCL+ 01
112 RCL 00
113 Y^X
114 2
115 ÷
116 +
117 1
118 RCL 00
119 +/-
120 STO 00
121 -
122 2
123 LN
124 ×
125 E^X-1
126 +/-
127 ÷
128 4
129 RCL- 00
130 44E-13
131 ×
132 X<0?
133 CLX
134 -
135 GTO 99
136>LBL 97
137 ENTER
138 ENTER
139 ENTER
140 -5.872675E-6
141 ×
142 9.740462E-5
143 +
144 ×
145 3.4213951E-4
146 -
147 ×
148 0.00484515533
149 -
150 ×
151 0.07281584734
152 +
153 ×
154 7.21566494432E-3
155 +
156>LBL 98
157 0.57
158 +
159 X<>Y
160 1/X
161 +
162>LBL 99
163 RCL 00
164 X<>Y
165 .END.

At least the results look quite good. Since x is returned in the upper display line you can see if a value very close to zero has been replaced by 1E–20.

Dieter

[Gerson W. Barbosa]

(08-04-2017 07:26 PM)Dieter Wrote:  [quote='Gerson W. Barbosa' pid='76831' dateline='1501861065']
Back to the HP-41:

Code:

 01▸LBL "ZETA"
 02 X>0?
 03 GTO 00
 04 X=0?
 05 GTO 00
 ...

Just one remark: Testing condition1 OR condition2 can be done with a simple trick from the olden days. Simply have the inverse of condition1 followed by condition2.

Code:

 01▸LBL "ZETA"
 02 X≠0?
 03 X>0?
 04 GTO 00
 ...

Code:


 02 X<=0?
 03 X=0?
 04 GTO OO
 ...

Lesson learned, thanks! I knew there was a better way, but I wouldn't remember. In the olden days I had no 41, but I did have a 15C which had all the 12 conditional tests, Gamma, complex mode and a few other thing the 41 lacked :-)

(08-04-2017 07:26 PM)Dieter Wrote:  Edit: here is the current version. It addresses the x~0 problem in a ...err... "creative" way in that the smallest |x| is set to 1E–20. #-)

Please note that this version is intended for high-precision environments like Free42 and may not perform as well on a regular 42s.

Code:

 00 { 416-Byte Prgm }
 01>LBL "ZETA"
 02 0.5
 03 X<>Y
 04 X>=Y?
 05 GTO 77
 06 ABS
 07 1E-20
 08 X<Y?
 09 X<>Y
 10 LASTX
 11 SIGN
 12 ×
 13 1
 14 X<>Y
 15 -
 16 XEQ 77
 17 PI
 18 STO+ ST X
 19 RCL 00
 20 Y^X
 21 ÷
 22 RCL 00
 23 GAMMA
 24 ×
 25 1
 26 ASIN
 27 RCL× 00
 28 COS
 29 ×
 30 STO+ ST X
 31 1
 32 RCL- 00
 33 STO 00
 34 Rv
 35 GTO 99
 36>LBL 77
 37 STO 00
 38 -1
 39 RCL+ 00
 40 1/X
 41 LASTX
 42 X<0?
 43 GTO 97
 44 2
 45 RCL 00
 46 X>Y?
 47 GTO 96
 48 LASTX
 49 LASTX
 50 LASTX
 51 -1.276E-8
 52 ×
 53 7.05133E-6
 54 -
 55 ×
 56 9.721157E-5
 57 +
 58 ×
 59 3.4243368E-4
 60 -
 61 ×
 62 0.00484515482
 63 -
 64 ×
 65 0.07281584288
 66 +
 67 ×
 68 7.215664988E-3
 69 +
 70 GTO 98
 71>LBL 96
 72 100
 73 RCL 00
 74 SQRT
 75 RCL× 00
 76 ÷
 77 5
 78 +
 79 IP
 80 STO+ ST X
 81 STO 01
 82 RCL 00
 83 +/-
 84 STO 00
 85 CLX
 86>LBL 95
 87 RCL ST Y
 88 RCL 00
 89 Y^X
 90 -
 91 +/-
 92 DSE ST Y
 93 GTO 95
 94 RCL 00
 95 STO+ ST X
 96 1
 97 -
 98 RCL 01
 99 X^2
100 24
101 ×
102 ÷
103 1
104 RCL- 00
105 8
106 ÷
107 RCL÷ 01
108 +
109 0.5
110 +
111 RCL+ 01
112 RCL 00
113 Y^X
114 2
115 ÷
116 +
117 1
118 RCL 00
119 +/-
120 STO 00
121 -
122 2
123 LN
124 ×
125 E^X-1
126 +/-
127 ÷
128 4
129 RCL- 00
130 44E-13
131 ×
132 X<0?
133 CLX
134 -
135 GTO 99
136>LBL 97
137 ENTER
138 ENTER
139 ENTER
140 -5.872675E-6
141 ×
142 9.740462E-5
143 +
144 ×
145 3.4213951E-4
146 -
147 ×
148 0.00484515533
149 -
150 ×
151 0.07281584734
152 +
153 ×
154 7.21566494432E-3
155 +
156>LBL 98
157 0.57
158 +
159 X<>Y
160 1/X
161 +
162>LBL 99
163 RCL 00
164 X<>Y
165 .END.

At least the results look quite good. Since x is returned in the upper display line you can see if a value very close to zero has been replaced by 1E–20.

This has replaced what I had in Free42 (per Natural Selection).

Gerson.

[Dieter]

(08-05-2017 04:36 AM)Gerson W. Barbosa Wrote:  

Code:

 02 X<=0?
 03 X=0?
 04 GTO OO
 ...

I see you get the idea. ;-)

(08-05-2017 04:36 AM)Gerson W. Barbosa Wrote:  Lesson learned, thanks! I knew there was a better way, but I wouldn't remember.

Starting quite exactly five years ago, this thread in the old forum may refresh your memories in this regard.

(08-05-2017 04:36 AM)Gerson W. Barbosa Wrote:  In the olden days I had no 41, but I did have a 15C which had all the 12 conditional tests, Gamma, complex mode and a few other thing the 41 lacked :-)

Hmmm... then I think of a complete set of flag tests as provided by the '41. ;-) Here the old tricks may help for the one or other elegant solution in a 15C program. On the 67/97 I love things like F2? F2? for a FC?C test, or three consecutive tests for more obscure applications. Or even four, like in the linked thread. I think I remember there also was an artice in an HP Journal issue...

(08-05-2017 04:36 AM)Gerson W. Barbosa Wrote:  This has replaced what I had in Free42 (per Natural Selection).

Fine. Let me know if you find any cases where – at least for x>0,5 – the 12th digit is off by more than 1 unit. I hope you don't. ;-)
With the applied correction term (line 128 ff.) the error for x>2 should be less than 1 unit in the 13th place.

Dieter