RE: Riemann's Zeta Function - another approach (RPL) - Dieter - 08-07-2017 08:52 AM
(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. #-)
There is a better solution.
Since \(\zeta '(0) = -\ln\sqrt{2 \pi}\), for x close to zero \(\zeta(x)\ \approx -(0.5+x\cdot\ln\sqrt{2 \pi})\). So the program may simply switch to this as soon as |x| is sufficiently small. For the given program a good threshold is at 2 E–11. So if |x| is larger than this, use the regular method, and all |x| up to this limit may use the simple formula above. This is what I now have implemented in my Free42 program.
Here is the modified listing. The changes are in the first lines before the XEQ 77 call, and the STO 00 after LBL 77 has been omitted.
Edit: changed the threshold in line 08 from 1E–11 to 2E–11.
Code:
00 { 434-Byte Prgm }
01>LBL "ZETA"
02 STO 00
03 0.5
04 X<>Y
05 X>=Y?
06 GTO 77
07 ABS
08 2E-11
09 X<Y?
10 GTO 00
11 PI
12 STO+ ST X
13 SQRT
14 LN
15 RCL× 00
16 0.5
17 +
18 +/-
19 GTO 99
20>LBL 00
21 1
22 RCL- 00
23 STO 00
24 XEQ 77
25 PI
26 STO+ ST X
27 RCL 00
28 Y^X
29 ÷
30 RCL 00
31 GAMMA
32 ×
33 1
34 ASIN
35 RCL× 00
36 COS
37 ×
38 STO+ ST X
39 1
40 RCL- 00
41 STO 00
42 Rv
43 GTO 99
44>LBL 77
45 -1
46 RCL+ 00
47 1/X
48 LASTX
49 X<0?
50 GTO 97
51 2
52 RCL 00
53 X>Y?
54 GTO 96
55 LASTX
56 LASTX
57 LASTX
58 -1.276E-8
59 ×
60 7.05133E-6
61 -
62 ×
63 9.721157E-5
64 +
65 ×
66 3.4243368E-4
67 -
68 ×
69 0.00484515482
70 -
71 ×
72 0.07281584288
73 +
74 ×
75 7.215664988E-3
76 +
77 GTO 98
78>LBL 96
79 100
80 RCL 00
81 SQRT
82 RCL× 00
83 ÷
84 5
85 +
86 IP
87 STO+ ST X
88 STO 01
89 RCL 00
90 +/-
91 STO 00
92 CLX
93>LBL 95
94 RCL ST Y
95 RCL 00
96 Y^X
97 -
98 +/-
99 DSE ST Y
100 GTO 95
101 RCL 00
102 STO+ ST X
103 1
104 -
105 RCL 01
106 X^2
107 24
108 ×
109 ÷
110 1
111 RCL- 00
112 8
113 ÷
114 RCL÷ 01
115 +
116 0.5
117 +
118 RCL+ 01
119 RCL 00
120 Y^X
121 2
122 ÷
123 +
124 1
125 RCL 00
126 +/-
127 STO 00
128 -
129 2
130 LN
131 ×
132 E^X-1
133 +/-
134 ÷
135 4
136 RCL- 00
137 44E-13
138 ×
139 X<0?
140 CLX
141 -
142 GTO 99
143>LBL 97
144 ENTER
145 ENTER
146 ENTER
147 -5.872675E-6
148 ×
149 9.740462E-5
150 +
151 ×
152 3.4213951E-4
153 -
154 ×
155 0.00484515533
156 -
157 ×
158 0.07281584734
159 +
160 ×
161 7.21566494432E-3
162 +
163>LBL 98
164 0.57
165 +
166 X<>Y
167 1/X
168 +
169>LBL 99
170 RCL 00
171 X<>Y
172 END
And again: please note that this version is intended for high-precision environments like Free42 and may not perform well on a regular 42s.
Dieter
RE: Riemann's Zeta Function - another approach (RPL) - Gerson W. Barbosa - 08-08-2017 04:05 AM
(08-07-2017 08:52 AM)Dieter Wrote: (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. #-)
There is a better solution.
Since \(\zeta '(0) = -\ln\sqrt{2 \pi}\), for x close to zero \(\zeta(x)\ \approx -(0.5+x\cdot\ln\sqrt{2 \pi})\). So the program may simply switch to this as soon as |x| is sufficiently small. For the given program a good threshold is at about 1E–11. So if |x| is larger than this, use the regular method, and all |x| up to this limit may use the simple formula above. This is what I now have implemented in my Free42 program.
Here is the modified listing. The changes are in the first lines before the XEQ 77 call, and the STO 00 after LBL 77 has been omitted.
Code:
00 { 434-Byte Prgm }
01>LBL "ZETA"
02 STO 00
03 0.5
04 X<>Y
05 X>=Y?
06 GTO 77
07 ABS
08 1E-11
09 X<Y?
10 GTO 00
11 PI
12 STO+ ST X
13 SQRT
14 LN
15 RCL× 00
16 0.5
17 +
18 +/-
19 GTO 99
20>LBL 00
21 1
22 RCL- 00
23 STO 00
24 XEQ 77
25 PI
26 STO+ ST X
27 RCL 00
28 Y^X
29 ÷
30 RCL 00
31 GAMMA
32 ×
33 1
34 ASIN
35 RCL× 00
36 COS
37 ×
38 STO+ ST X
39 1
40 RCL- 00
41 STO 00
42 Rv
43 GTO 99
44>LBL 77
45 -1
46 RCL+ 00
47 1/X
48 LASTX
49 X<0?
50 GTO 97
51 2
52 RCL 00
53 X>Y?
54 GTO 96
55 LASTX
56 LASTX
57 LASTX
58 -1.276E-8
59 ×
60 7.05133E-6
61 -
62 ×
63 9.721157E-5
64 +
65 ×
66 3.4243368E-4
67 -
68 ×
69 0.00484515482
70 -
71 ×
72 0.07281584288
73 +
74 ×
75 7.215664988E-3
76 +
77 GTO 98
78>LBL 96
79 100
80 RCL 00
81 SQRT
82 RCL× 00
83 ÷
84 5
85 +
86 IP
87 STO+ ST X
88 STO 01
89 RCL 00
90 +/-
91 STO 00
92 CLX
93>LBL 95
94 RCL ST Y
95 RCL 00
96 Y^X
97 -
98 +/-
99 DSE ST Y
100 GTO 95
101 RCL 00
102 STO+ ST X
103 1
104 -
105 RCL 01
106 X^2
107 24
108 ×
109 ÷
110 1
111 RCL- 00
112 8
113 ÷
114 RCL÷ 01
115 +
116 0.5
117 +
118 RCL+ 01
119 RCL 00
120 Y^X
121 2
122 ÷
123 +
124 1
125 RCL 00
126 +/-
127 STO 00
128 -
129 2
130 LN
131 ×
132 E^X-1
133 +/-
134 ÷
135 4
136 RCL- 00
137 44E-13
138 ×
139 X<0?
140 CLX
141 -
142 GTO 99
143>LBL 97
144 ENTER
145 ENTER
146 ENTER
147 -5.872675E-6
148 ×
149 9.740462E-5
150 +
151 ×
152 3.4213951E-4
153 -
154 ×
155 0.00484515533
156 -
157 ×
158 0.07281584734
159 +
160 ×
161 7.21566494432E-3
162 +
163>LBL 98
164 0.57
165 +
166 X<>Y
167 1/X
168 +
169>LBL 99
170 RCL 00
171 X<>Y
172 END
Very nice!
I would only make a minor modification. Not that one less square root evaluation really might make a difference in execution time in Free42:
Code:
…
13 LN
14 RCL× 00
15 RCL× ST T
16 RCL+ ST T
17 +/-
…
Gerson.
RE: Riemann's Zeta Function - another approach (RPL) - Dieter - 08-08-2017 06:46 PM
(08-08-2017 04:05 AM)Gerson W. Barbosa Wrote: Very nice!
Wait, it even gets a tiny bit nicer. ;-)
(08-08-2017 04:05 AM)Gerson W. Barbosa Wrote: I would only make a minor modification. Not that one less square root evaluation really might make a difference in execution time in Free42:
Code:
…
13 LN
14 RCL× 00
15 RCL× ST T
16 RCL+ ST T
17 +/-
…
Ah, great. It's not beause of the sqrt, but I love this elegant stack content reuse of the 0.5 that still sits in T. ;-)
But there even is one more tweak. The second term of the Zeta Taylor series is quite exactly –x², i.e. \(\zeta(x)\ \approx -(0.5+x^2+x\cdot\ln\sqrt{2 \pi})\).
And this can be implemented with just one more step:
Code:
…
13 LN
14 RCL× ST T
15 RCL+ 00
16 RCL× 00
17 RCL+ ST T
18 +/-
…
This way the threshold in line 08 should be somewhere near 1 E–8. The perfect value differs slightly for negative resp. positive x, but 8 E–9 works very nicely for both.
BTW, if you want to stick to the original implementation the perfect switchpoint seems to be close to 2 E–11, so this value should be used there. I edited my original post accordingly. But we're talking about differences beyond the 20th digit here. ;-)
Dieter
RE: Riemann's Zeta Function - another approach (RPL) - Gerson W. Barbosa - 08-12-2017 07:03 PM
(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
...
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
Just for the record, here is the new HP-41 listing with your suggestions:
Code:
01▸LBL "ZETA"
02 X<=0?
03 X=0?
04 GTO 00
05 CHS
06 1
07 +
08 STO 02
09 XEQ "GAMMA"
10 STO 03
11 RCL 02
12 XEQ 00
13 RCL 03
14 *
15 PI
16 ST+ X
17 RCL 02
18 Y^X
19 /
20 1
21 ASIN
22 RCL 02
23 *
24 COS
25 *
26 ST+ X
27 GTO 99
28▸LBL 00
29 STO 00
30 1
31 -
32 1/X
33 LASTX
34 X<0?
35 GTO 97
36 2
37 RCL 00
38 X>Y?
49 GTO 96
40 LASTX
41 LASTX
42 LASTX
43 -1.276E -8
44 *
45 7.05133E -6
46 -
47 *
48 9.721157E -5
49 +
50 *
51 3.4243368E -4
52 -
53 *
54 4.84515482E -3
55 -
56 *
57 7.281584288E -2
58 +
59 *
60 7.215664988E -3
61 +
62 GTO 98
63▸LBL 96
64 24
65 RCL 00
66 /
67 4
68 +
69 INT
70 ST+ X
71 STO 01
72 RCL 00
73 CHS
74 STO 00
75 CLX
76▸LBL 01
77 RCL Y
78 RCL 00
79 Y^X
80 -
81 CHS
82 DSE Y
83 GTO 01
84 RCL 00
85 ST+ X
86 1
87 -
88 RCL 01
89 X^2
90 24
91 *
92 /
93 1
94 RCL 00
95 -
96 8
97 /
98 RCL 01
99 /
100 +
101 .5
102 +
103 RCL 01
104 +
105 RCL 00
106 Y^X
107 2
108 /
109 +
110 1
111 RCL 00
112 +
113 2
114 LN
115 *
116 E^X-1
117 CHS
118 /
119 GTO 99
120▸LBL 97
121 ENTER^
122 ENTER^
123 ENTER^
124 -8.4715E -7
125 *
126 7.51334E -6
127 -
128 *
129 9.609657E -5
130 +
131 *
132 3.42683396E -4
133 -
134 *
135 4.84527616E -3
136 -
137 *
138 7.281583446E -2
139 +
140 *
141 7.215664464E -3
142 +
143▸LBL 98
144 RDN
145 1/X
146 INT
147 ST* Z
148 SIGN
149 STO/ X
150 STO- Z
151 X<> L
152 RDN
153 /
154 -
155 R^
156 .57
157 +
158 +
159▸LBL 99
160 END
The primitive Gamma implementation I've been using is slightly faster:
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 DSE 00
18 ABS
19 DSE Y
20 GTO 00
21 STO 00
22▸LBL 01
23 LASTX
24 ENTER^
25 ENTER^
26 ENTER^
27 1.604589926 E-2
28 *
29 2.641400819 E-1
30 -
31 *
32 1.96580515
33 +
34 *
35 8.729327662
36 -
37 *
38 25.69590147
39 +
40 *
41 52.63472435
42 -
43 *
44 76.53826433
45 +
46 *
47 78.92639573
48 -
49 *
50 56.50761084
51 +
52 *
53 26.37243835
54 -
55 *
56 7.203398486
57 +
58 RCL 00
59 *
60 END
And here is the new table:
41 XEQ ZETA --> 1.000000000 ( 8.24 s)
25 R/S --> 1.000000030 ( 8.28 s)
3 R/S --> 1.202056903 (17.89 s)
2.001 R/S --> 1.643997513 (2) (22.33 s)
2 R/S --> 1.644934067 ( 4.04 s)
1.5 R/S --> 2.612375349 ( 4.08 s)
0.5 R/S --> -1.460354509 ( 4.74 s)
0 R/S --> -0.500000000 ( 4.03 s)
-0.5 R/S --> -0.2078862450 (250) ( 9.72 s)
-1 R/S --> -0.08333333384 (33) (10.07 s)
-1.001 R/S --> -0.08316803696 (723) (28.32 s)
-1.5 R/S --> -0.02548520436 (190) (25.38 s)
-2 R/S --> 0.00000000000 (24.84 s)
-3 R/S --> 0.008333333384 (33) (22.09 s)
-5 R/S --> -0.003968253990 (68) (20.24 s)
-7 R/S --> 0.004166666686 (67) (19.39 s)
-15.16 R/S --> 0.4964873534 (85) (18.99 s)
-33.34 R/S --> -1.924684098E10 (152) (22.15 s)
-41.42 R/S --> -3.506595602E16 (584) (24.03 s)
-48.49 R/S --> -3.653091058E22 (22) (25.98 s)
-58.59 R/S --> 1.136304829E32 (792) (28.54 s)
-67.97 R/S --> 1.832460467E40 (1182) (31.02 s)
Times obtained from the HP-41CX built-in stopwatch
RE: Riemann's Zeta Function - another approach (RPL) - Gerson W. Barbosa - 05-08-2020 04:41 PM
Same as before, except for the shorter (90 bytes, the other verstion was 191 bytes long and a bit inaccurate near zero) and more accurate Gamma funcion.
To Dieter, who surely would improve this "improved" Gamma program as he did to Zeta, were he still among us. This nice Zeta implementation is mostly his accomplishment.
HP-41CV/CX
ζ(x)
Code:
01 LBL "ZETA"
02 X<=0?
03 X=0?
04 GTO 00
05 CHS
06 1
07 +
08 STO 02
09 XEQ "GAMMA"
10 STO 03
11 RCL 02
12 XEQ 00
13 RCL 03
14 *
15 PI
16 ST+ X
17 RCL 02
18 Y^X
19 /
20 1
21 ASIN
22 RCL 02
23 *
24 COS
25 *
26 ST+ X
27 GTO 99
28 LBL 00
29 STO 00
30 1
31 -
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 E-8
44 *
45 7.05133 E-6
46 -
47 *
48 9.721157 E-5
49 +
50 *
51 3.4243368 E-4
52 -
53 *
54 4.84515482 E-3
55 -
56 *
57 7.281584288 E-2
58 +
59 *
60 7.215664988 E-3
61 +
62 GTO 98
63 LBL 96
64 24
65 RCL 00
66 /
67 4
68 +
69 INT
70 ST+ X
71 STO 01
72 RCL 00
73 CHS
74 STO 00
75 CLX
76 LBL 01
77 RCL Y
78 RCL 00
79 Y^X
80 -
81 CHS
82 DSE Y
83 GTO 01
84 RCL 00
85 ST+ X
86 1
87 -
88 RCL 01
89 X^2
90 24
91 *
92 /
93 1
94 RCL 00
95 -
96 8
97 /
98 RCL 01
99 /
100 +
101 .5
102 +
103 RCL 01
104 +
105 RCL 00
106 Y^X
107 2
108 /
109 +
110 1
111 RCL 00
112 +
113 2
114 LN
115 *
116 E^X-1
117 CHS
118 /
119 GTO 99
120 LBL 97
121 ENTER
122 ENTER
123 ENTER
124 -8.4715 E-7
125 *
126 7.51334 E-6
127 -
128 *
129 9.609657 E-5
130 +
131 *
132 3.42683396 E-4
133 -
134 *
135 4.84527616 E-3
136 -
137 *
138 7.281583446 E-2
139 +
140 *
141 7.215664464 E-3
142 +
143 LBL 98
144 RDN
145 1/X
146 INT
147 ST* Z
148 SIGN
149 ST/ X
150 ST- Z
151 X<> L
152 RDN
153 /
154 -
155 R^
156 .57
157 +
158 +
159 LBL 99
160 END
362 BYTES
Γ(x), x ≥ 0
Code:
01 LBL "GAMMA"
02 1
03 STO 00
04 -
05 4
06 STO Z
07 X<=Y?
08 GTO 03
09 +
10 RCL X
11 ISG X
12 LBL 02
13 1
14 -
15 ST* 00
16 DSE Z
17 GTO 02
18 LBL 03
19 X<>Y
20 4.13333
21 2.215
22 RCL Z
23 *
24 .4875
25 -
26 R^
27 X^2
28 /
29 +
30 12
31 /
32 +
33 72
34 *
35 1/X
36 6
37 1/X
38 +
39 +
40 360
41 D-R
42 *
43 SQRT
44 X<>Y
45 1
46 E^X
47 /
48 R^
49 Y^X
50 *
51 RCL 00
52 /
53 END
90 BYTES
41 XEQ ZETA --> 1.000000000 ( 8.3 s)
25 R/S --> 1.000000030 ( 8.3 s)
3 R/S --> 1.202056903 (17.5 s)
2.001 R/S --> 1.643997513 (2) (21.7 s)
2 R/S --> 1.644934067 ( 5.0 s)
1.5 R/S --> 2.612375349 ( 4.3 s)
0.5 R/S --> -1.460354509 ( 4.3 s)
0 R/S --> -0.500000000 ( 4.3 s)
-0.5 R/S --> -0.2078862256 (0) (10.3 s)
-1 R/S --> -0.08333333344 (33) ( 9.8 s)
-1.001 R/S --> -0.08316803746 (46) (27.4 s)
-1.5 R/S --> -0.02548520190 (24.6 s)
-2 R/S --> 0.00000000000 (23.0 s)
-3 R/S --> 0.008333333350 (33) (20.9 s)
-5 R/S --> -0.003968253966 (8) (17.7 s)
-7 R/S --> 0.004166666668 (7) (16.4 s)
-15.16 R/S --> 0.4964873582 (2) (14.5 s)
-33.34 R/S --> -1.924684152E10 (13.1 s)
-41.42 R/S --> -3.506595630E16 (584) (13.1 s)
-48.49 R/S --> -3.653091072E22 (22) (13.2 s)
-58.59 R/S --> 1.136304789E32 (92) (13.2 s)
-67.97 R/S --> 1.832461183E40 (2) (13.3 s)
Times on my HP-41CV
————
P.S.:
For a full-range Γ(x+1) implementation please refer to
Γ(x+1) [HP-41C].
————
RE: Riemann's Zeta Function - another approach (RPL) - rprosperi - 05-08-2020 04:47 PM
(05-08-2020 04:41 PM)Gerson W. Barbosa Wrote: ...To Dieter, who surely would improve this "improved" Gamma program as he did to Zeta, were he still among us. This nice Zeta implementation is mostly his accomplishment.
Nice callout Gerson. I think of Dieter every time I see a post of someone improve and shorten their own code...
I guess he knows all the tricks now...
RE: Riemann's Zeta Function - another approach (RPL) - Paul Dale - 05-09-2020 12:11 AM
Agreed, a nice eulogy. I miss Dieter's insightful and erudite input.
Pauli
|