Riemann's Zeta Function - another approach (RPL)

[Dieter]

(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

[Gerson W. Barbosa]

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

[Dieter]

(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

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

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

[Gerson W. Barbosa]

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

————

[rprosperi]

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

[Paul Dale]

Agreed, a nice eulogy. I miss Dieter's insightful and erudite input.

Pauli