Zeta

Riemann Zeta Function for the HP-41

Overview

1°) Real Arguments

a) Elementary Method
b) Borwein Algorithm

2°) Complex Arguments
3°) Quaternion Arguments

a) Elementary Method
b) Borwein Algorithm

1°) Real Arguments

a) Elementary Method

-The Riemann Zeta function is defined by the series: Zeta(x) = 1 + 1/2^x + 1/3^x + ........ + 1/n^x + ...... ( for x > 1)
but the following program uses the formula: Zeta(x) = ( 1 - 1/2^x + 1/3^x - 1/4^x + ...... ) / ( 1 - 2^1-x ) where roundoff errors are smaller.

-If x < 0.5 , we have: Zeta(x) = GAM(1-x) 2^x PI^x-1 Sin(PI.x/2) Zeta(1-x)

Data Registers: R00 thru R03 ( temp )
Flags: /
Subroutine: "GAM" ( cf "Gamma Function for the HP-41" ) if x < 1 only.

01 LBL "ZETA"
02 1
03 X<>Y
04 X>Y?
05 GTO 00
06 STO 02
07 -
08 STO 03
09 XEQ 00
10 1
11 ASIN
12 RCL 02

13 *
14 SIN
15 *
16 PI
17 RCL 03
18 Y^X
19 /
20 2
21 RCL 02
22 Y^X
23 *
24 X<> 03

25 XEQ "GAM"
26 RCL 03
27 *
28 RTN
29 LBL 00
30 CHS
31 STO 01
32 1
33 STO 00
34 ENTER^
35 LBL 01
36 CLX

37 SIGN
38 RCL 00
39 +
40 STO 00
41 RCL 01
42 Y^X
43 -
44 CHS
45 LASTX
46 RND
47 X#0?
48 GTO 01

49 SIGN
50 2
51 RCL 01
52 Y^X
53 ST+ X
54 -
55 /
56 ABS
57 END

( 76 bytes / SIZE 004 )

STACK	INPUT	OUTPUT
X	x	Zeta(x)

Examples:

FIX 7 3 XEQ "ZETA" returns Zeta(3) = 1.2020569 ( in 3m14s )
FIX 9 3 R/S ------- Zeta(3) = 1.202056903 ( in 15m12s )

FIX 9 -7.49 R/S ------- Zeta(-7.49) = 0.003312040169 ( in 15s )

Notes:

-The accuracy is determined by the display format ( line 46 ).
-Execution time and Zeta(x) tend to infinity as x tends to 1.

b) Borwein Algorithm

-This second version uses the method described in reference [3]

Data Registers: R00 to R07: temp
Flags: /
Subroutine: "1/G+" ( cf "Gamma function for the HP-41" )

-Line 69 is a synthetic TEXT0 instruction ( NOP )

01 LBL "ZETA"
02 X#0?
03 GTO 00
04 .5
05 CHS
06 RTN
07 LBL 00
08 STO 01
09 STO 07
10 .5
11 X<=Y?
12 GTO 00
13 SIGN
14 X<>Y
15 -
16 STO 01
17 XEQ 00
18 1

19 ASIN
20 RCL 07
21 *
22 SIN
23 *
24 PI
25 ST/ Y
26 ST+ X
27 RCL 07
28 Y^X
29 *
30 X<> 01
31 XEQ "1/G+"
32 ST/ 01
33 X<> 01
34 RTN
35 LBL 00
36 14

37 STO 00
38 STO 03
39 SIGN
40 STO 02
41 STO 04
42 STO 05
43 CHS
44 STO 06
45 CLX
46 LBL 01
47 RCL 05
48 RCL 00
49 RCL 01
50 CHS
51 Y^X
52 *
53 RCL 06
54 CHS

55 STO 06
56 *
57 +
58 RCL 00
59 ENTER^
60 ST+ Y
61 ST* Y
62 -
63 RCL 04
64 *
65 RCL 03
66 X^2
67 RCL 00
68 DSE X
69 ""
70 X^2
71 -
72 ST+ X

73 /
74 STO 04
75 ST+ 05
76 X<>Y
77 DSE 00
78 GTO 01
79 RCL 05
80 /
81 2
82 LN
83 1
84 RCL 01
85 -
86 *
87 E^X-1
88 /
89 END

( 122 bytes / SIZE 008 )

STACK	INPUT	OUTPUT
X	x	Zeta(x)

Examples:

      3    XEQ "ZETA"   >>>>      Zeta(3)    = 1.202056903                   ---Execution time = 21s---
-7.49 XEQ "ZETA"   >>>>    Zeta(-7.49) = 0.003312040169            ---Execution time = 24s---
    1.1   XEQ "ZETA" >>>>      Zeta(1.1) = 10.58444847                   ---Execution time = 21s---

Notes:

Zeta(0) = -1/2
-Borwein Algorithm is really fantastic !
-You get Zeta(1.001) in 21 seconds, which seemed impossible before...

2°) Complex Arguments

-"ZETAZ" also uses Borwein algorithm if Re(z) > 0.5

-If Re(z) < 0.5 , we've used Zeta(z) = Zeta(1-z) Pi^z-1/2 Gamma((1-z)/2) / Gamma(z/2)

-The other formula may also be used but it seems to produce more roundoff-errors.

Data Registers: R00 to R11: temp      When the program stops, the result is in R04-R05
Flags: /
Subroutines:   "GAMZ+" ( cf "Gamma Function for the HP-41" )
                           "Z*Z" "Z/Z" "Z^Z" ( cf "Complex Functions for the HP-41" )

-Line 128 is a synthetic TEXT0 instruction ( NOP )

01 LBL "ZETAZ"
02 X=0?
03 X#Y?
04 GTO 00
05 .5
06 CHS
07 RTN
08 LBL 00
09 X<>Y
10 STO 07
11 STO 11
12 X<>Y
13 STO 06
14 STO 10
15 .5
16 X<=Y?
17 GTO 00
18 SIGN
19 X<>Y
20 -
21 STO 06
22 RCL 07
23 CHS
24 STO 07
25 XEQ 00
26 RCL 11
27 CHS
28 1
29 RCL 10
30 -

31 2
32 ST/ Z
33 /
34 XEQ "GAMZ+"
35 RCL 05
36 RCL 04
37 XEQ "Z*Z"
38 STO 04
39 X<>Y
40 STO 05
41 PI
42 RCL 11
43 RCL 10
44 .5
45 ST* 10
46 ST* 11
47 -
48 0
49 RDN
50 XEQ "Z^Z"
51 RCL 05
52 RCL 04
53 XEQ "Z*Z"
54 STO 04
55 X<>Y
56 STO 05
57 RCL 11
58 RCL 10
59 XEQ "GAMZ+"
60 RCL 05

61 RCL 04
62 GTO 02
63 LBL 00
64 PI
65 2
66 /
67 RCL 07
68 ABS
69 ST* Y
70 ST+ X
71 LN1+X
72 +
73 3 E10
74 LN
75 +
76 8
77 SQRT
78 3
79 +
80 LN
81 /
82 INT
83 1
84 +
85 STO 00
86 STO 02
87 LASTX
88 STO 01
89 STO 03
90 STO 04

91 CHS
92 X<>Y
93 Y^X
94 CHS
95 STO 05
96 CLX
97 STO 08
98 STO 09
99 LBL 01
100 CLST
101 RCL 00
102 RCL 07
103 CHS
104 RCL 06
105 CHS
106 XEQ "Z^Z"
107 RCL 04
108 RCL 05
109 CHS
110 STO 05
111 *
112 ST* Z
113 *
114 ST+ 08
115 X<>Y
116 ST+ 09
117 RCL 00
118 ENTER^
119 ST+ Y
120 ST* Y

121 -
122 RCL 03
123 *
124 RCL 02
125 X^2
126 RCL 00
127 DSE X
128 ""
129 X^2
130 -
131 ST+ X
132 /
133 STO 03
134 ST+ 04
135 DSE 00
136 GTO 01
137 RCL 04
138 ST/ 08
139 ST/ 09
140 2
141 LN
142 1
143 RCL 06
144 -
145 *
146 E^X-1
147 STO 01
148 RCL 07
149 CHS
150 2

151 LN
152 *
153 1
154 RAD
155 P-R
156 ENTER^
157 DEG
158 1
159 -
160 RCL 01
161 ST* Z
162 X<> T
163 ST* T
164 ST+ T
165 RDN
166 +
167 RCL 09
168 RCL 08
169 LBL 02
170 R^
171 R^
172 XEQ "Z/Z"
173 STO 04
174 X<>Y
175 STO 05
176 X<>Y
177 END

( 251 bytes / SIZE 012 )

STACK	INPUTS	OUTPUTS
Y	y	y'
X	x	x'

Zeta(x+i.y) = x' + i.y'

Example:

7 ENTER^
-6 XEQ "ZETAZ" >>>> -4.135010534
RDN -1.198478879

-So, Zeta(-6+7.i) = -4.135010534 - 1.198478879 i

3°) Quaternion Arguments

a) Elementary Method

-"ZETAQ" uses Zeta(q) = 1 + 1/2^q + 1/3^q + ........ + 1/n^q + ...... ( for Re(q) > 1)

-And if Re(q) < 0 , Zeta(q) = Zeta(1-q) Pi^q-1/2 Gamma((1-q)/2) / Gamma(q/2)

Data Registers: R00 to R24: temp
Flags: /
Subroutines: "GAMQ" ( cf "Quaternionic Special Functions for the HP-41" )
"Q*Q" "1/Q" "E^Q" ( cf "Quaternions for the HP-41" )

-Line 16 is a three-byte GTO 00

01 LBL "ZETAQ"
02 STO 05
03 STO 17
04 RDN
05 STO 06
06 STO 18
07 RDN
08 STO 07
09 STO 19
10 X<>Y
11 STO 08
12 STO 20
13 RCL 05
14 1
15 X<Y?
16 GTO 00
17 CHS
18 ST* 05
19 ST* 06
20 ST* 07
21 ST* 08
22 ST- 05
23 XEQ 00
24 STO 21
25 RDN
26 STO 22
27 RDN
28 STO 23
29 X<>Y
30 STO 24

31 2
32 ST/ 05
33 ST/ 06
34 ST/ 07
35 ST/ 08
36 RCL 08
37 RCL 07
38 RCL 06
39 RCL 05
40 XEQ "GAMQ"
41 STO 05
42 RDN
43 STO 06
44 RDN
45 STO 07
46 X<>Y
47 STO 08
48 21.001004
49 REGMOVE
50 XEQ "Q*Q"
51 STO 21
52 RDN
53 STO 22
54 RDN
55 STO 23
56 X<>Y
57 STO 24
58 RCL 17
59 SIGN
60 RCL 20

61 RCL 19
62 RCL 18
63 2
64 ST/ L
65 ST/ Y
66 ST/ Z
67 ST/ T
68 X<> L
69 XEQ "GAMQ"
70 XEQ "1/Q"
71 STO 05
72 RDN
73 STO 06
74 RDN
75 STO 07
76 X<>Y
77 STO 08
78 21.001004
79 REGMOVE
80 XEQ "Q*Q"
81 STO 05
82 RDN
83 STO 06
84 RDN
85 STO 07
86 X<>Y
87 STO 08
88 .5
89 ST- 17
90 RCL 20

91 RCL 19
92 RCL 18
93 PI
94 LN
95 ST* 17
96 ST* T
97 ST* Z
98 *
99 RCL 17
100 XEQ "E^Q"
101 STO 01
102 RDN
103 STO 02
104 RDN
105 STO 03
106 X<>Y
107 STO 04
108 XEQ "Q*Q"
109 RTN
110 LBL 00
111 CLX
112 STO 02
113 STO 03
114 STO 04
115 SIGN
116 STO 00
117 STO 01
118 LBL 01
119 ISG 00
120 CLX

121 RCL 08
122 RCL 07
123 RCL 06
124 RCL 00
125 LN
126 CHS
127 SIGN
128 CLX
129 RCL 05
130 X<> L
131 ST* L
132 ST* Y
133 ST* Z
134 ST* T
135 X<> L
136 XEQ "E^Q"
137 STO 09
138 CLX
139 RCL 02
140 +
141 STO 02
142 LASTX
143 -
144 ABS
145 X<>Y
146 RCL 03
147 +
148 STO 03
149 LASTX
150 -

151 ABS
152 +
153 X<>Y
154 RCL 04
155 +
156 STO 04
157 LASTX
158 -
159 ABS
160 +
161 RCL 09
162 RCL 01
163 +
164 STO 01
165 LASTX
166 -
167 ABS
168 +
169 X#0?
170 GTO 01
171 RCL 04
172 RCL 03
173 RCL 02
174 RCL 01
175 END

( 288 bytes / SIZE 025 )

STACK	INPUTS	OUTPUTS
T	t	t'
Z	z	z'
Y	y	y'
X	x	x'

Zeta(x+i.y+j.z+k.t) = x' + i.y' + j.z' + k.t'

Example:

   7.9   ENTER^
   7.8   ENTER^
   7.7   ENTER^
-7.6   XEQ "ZETAQ"    >>>>     762.9852911                                    ---Execution time = 2mn---
                                      RDN      -14.26823437
                                      RDN      -14.45353627
                                      RDN      -14.63883820

-Thus, Zeta ( -7.6 + 7.7 i + 7.8 j + 7.9 k ) = 762.9852911 - 14.26823437 i - 14.45353627 j - 14.63883820 k

Notes:

-Lines 137 to 168 may be replaced by DSQ 1 ( cf "M-Code routines for hypercomplex numbers" )
-The execution time corresponds to this modification...

-Add RND after line 168 and the accuracy will be controlled by the display settings.

b) Borwein Algorithm

Data Registers: R00 to R28: temp
Flags: /
Subroutines: "GAMQ" ( cf "Quaternionic Special Functions for the HP-41" )
"Q*Q" "1/Q" "E^Q" ( cf "Quaternions for the HP-41" )

-Lines 30 and 137 are three-byte GTOs

01 LBL "ZETAQ"
02 STO 05
03 STO 17
04 X^2
05 X<>Y
06 STO 06
07 STO 18
08 X^2
09 +
10 X<>Y
11 STO 07
12 STO 19
13 X^2
14 +
15 X<>Y
16 STO 08
17 STO 20
18 X^2
19 +
20 X#0?
21 GTO 00
22 CLST
23 .5
24 CHS
25 RTN
26 LBL 00
27 RCL 05
28 .5
29 X<=Y?
30 GTO 00
31 SIGN
32 CHS
33 ST* 05
34 ST* 06
35 ST* 07
36 ST* 08
37 ST- 05
38 RCL 05
39 STO 25
40 RCL 06
41 STO 26
42 RCL 07
43 STO 27
44 RCL 08
45 STO 28
46 XEQ 00
47 STO 21
48 RDN

49 STO 22
50 RDN
51 STO 23
52 X<>Y
53 STO 24
54 RCL 28
55 RCL 27
56 RCL 26
57 2
58 ST/ 25
59 ST/ T
60 ST/ Z
61 /
62 RCL 25
63 XEQ "GAMQ"
64 STO 05
65 RDN
66 STO 06
67 RDN
68 STO 07
69 X<>Y
70 STO 08
71 RCL 21
72 STO 01
73 RCL 22
74 STO 02
75 RCL 23
76 STO 03
77 RCL 24
78 STO 04
79 XEQ "Q*Q"
80 STO 21
81 RDN
82 STO 22
83 RDN
84 STO 23
85 X<>Y
86 STO 24
87 RCL 20
88 RCL 19
89 RCL 18
90 RCL 17
91 SIGN
92 CLX
93 2
94 ST/ L
95 ST/ Y
96 ST/ Z

97 ST/ T
98 X<> L
99 XEQ "GAMQ"
100 XEQ "1/Q"
101 STO 05
102 RDN
103 STO 06
104 RDN
105 STO 07
106 X<>Y
107 STO 08
108 RCL 21
109 STO 01
110 RCL 22
111 STO 02
112 RCL 23
113 STO 03
114 RCL 24
115 STO 04
116 XEQ "Q*Q"
117 STO 01
118 RDN
119 STO 02
120 RDN
121 STO 03
122 X<>Y
123 STO 04
124 .5
125 ST- 17
126 PI
127 LN
128 ST* 17
129 ST* 18
130 ST* 19
131 ST* 20
132 RCL 20
133 RCL 19
134 RCL 18
135 RCL 17
136 XEQ "E^Q"
137 GTO 02
138 LBL 00
139 PI
140 2
141 /
142 RCL 06
143 X^2
144 RCL 07

145 X^2
146 RCL 08
147 X^2
148 +
149 +
150 SQRT
151 ST* Y
152 ST+ X
153 LN1+X
154 +
155 3 E10
156 LN
157 +
158 8
159 SQRT
160 3
161 +
162 LN
163 /
164 INT
165 1
166 +
167 STO 00
168 STO 02
169 LASTX
170 STO 01
171 STO 03
172 STO 04
173 CHS
174 X<>Y
175 Y^X
176 CHS
177 STO 13
178 CLX
179 STO 09
180 STO 10
181 STO 11
182 STO 12
183 LBL 01
184 RCL 00
185 LN
186 SIGN
187 RCL 08
188 CHS
189 RCL 07
190 CHS
191 RCL 06
192 CHS

193 RCL 05
194 CHS
195 X<> L
196 ST* L
197 ST* Y
198 ST* Z
199 ST* T
200 X<> L
201 XEQ "E^Q"
202 STO 14
203 RDN
204 STO 15
205 CLX
206 RCL 04
207 RCL 13
208 CHS
209 STO 13
210 *
211 ST* 14
212 ST* 15
213 ST* Z
214 *
215 ST+ 11
216 X<>Y
217 ST+ 12
218 RCL 14
219 ST+ 09
220 RCL 15
221 ST+ 10
222 RCL 00
223 ENTER^
224 ST+ Y
225 ST* Y
226 -
227 RCL 03
228 *
229 RCL 02
230 X^2
231 RCL 00
232 1
233 -
234 X^2
235 -
236 ST+ X
237 /
238 STO 03
239 ST+ 04
240 DSE 00

241 GTO 01
242 RCL 09
243 RCL 04
244 /
245 STO 01
246 RCL 10
247 LASTX
248 /
249 STO 02
250 RCL 11
251 LASTX
252 /
253 STO 03
254 RCL 12
255 LASTX
256 /
257 STO 04
258 RCL 08
259 RCL 07
260 RCL 06
261 2
262 LN
263 CHS
264 ST* 05
265 ST- 05
266 ST* T
267 ST* Z
268 *
269 RCL 05
270 XEQ "E^Q"
271 SIGN
272 ST* X
273 ST- L
274 X<> L
275 XEQ "1/Q"
276 LBL 02
277 STO 05
278 RDN
279 STO 06
280 RDN
281 STO 07
282 X<>Y
283 STO 08
284 XEQ "Q*Q"
285 END

( 432 bytes / SIZE 029 )

STACK	INPUTS	OUTPUTS
T	t	t'
Z	z	z'
Y	y	y'
X	x	x'

Zeta(x+i.y+j.z+k.t) = x' + i.y' + j.z' + k.t'

Examples:

   7.9   ENTER^
   7.8   ENTER^
   7.7   ENTER^
-7.6   XEQ "ZETAQ"    >>>>     762.9852917                          ---Execution time = 2m45s---
                                      RDN      -14.26823447
                                      RDN      -14.45353622
                                      RDN      -14.63883832

Zeta ( -7.6 + 7.7 i + 7.8 j + 7.9 k ) = 762.9852917 - 14.26823447 i - 14.45353622 j - 14.63883832 k

-Likewise, Zeta( 1.1 + 7 i + 8 j + 9 k ) = 0.389296030 - 0.027737071 i - 0.031699509 j - 0.035661948 k ( in 2m06s )

Notes:

-If q is very close to 1, there could be a loss of accuracy when 2^(1-q)-1 is executed:
-Replace lines 258 to 275 by

   1                   STO 09          RCL 16        ST* Y             RCL 10          CHS               ST* 06           RCL 06
   RCL 05         *                    *                   1                     E^X               RCL 16           ST* 07          RCL 05
   -                    STO 10          STO 09        -                     RCL 09          X=0?               ST* 08          DEG
   2                    E^X-1            RAD            +                    SIN                SIGN              RCL 08          XEQ "1/Q"
   LN                RCL 09          COS            STO 05           *                     /                     RCL 07

and add STO 16 after line 150.

-The last decimals may depend on the version of "E^Q" , "Q*Q" , .... that you are using.
-If Re(q) is large, it could be preferable to employ the 1st version of "ZETAQ".
-In all other cases, Borwein algorithm is unbeatable !

References:

[1] Abramowitz and Stegun , "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4
[2] http://functions.wolfram.com/
[3] P. Borwein - "An Efficient Algorithm for the Riemann Zeta Function"

hp41programs

Riemann Zeta Function for the HP-41