Riemann's Zeta Function (HP 50g) Message #1 Posted by Gerson W. Barbosa on 4 Dec 2012, 3:59 p.m.
This is a followup of a recent thread on the Zeta function . On the occasion Thomas Ritschel wrote an RPL program based on EulerMacLaurin summation, which was much faster and efficient than my previous attempts. The following is based on the same method, as described in the book Riemann's Zeta Function, by H. M. Edwards (pages 114 through 118).
Zeta:
%%HP: T(3)A(D)F(.);
\<< RCLF SWAP 22. SF RAD 1. CF DUP RE 0. <
IF
THEN 1. SF NEG 1. +
END 10. 8. 0. \> s n m b
\<< '\GS(k=1.,n1.,k^s)' EVAL n 1. s  ^ s 1.  / + n s NEG ^ 2. / + BNL 1. 2. m
FOR i GETI s i + 1.  GAMMA * n 1. i  s  ^ * s GAMMA / i ! / 'b' STO+ 2.
STEP b NIP NIP + 1. FS?
IF
THEN 1. s  DUP 2. / \pi * SIN 2. ROT ^ * * s DUP GAMMA \pi ROT NEG ^ * *
END
\>> SWAP STOF
\>>
BNL:
%%HP: T(3)A(D)F(.);
{ .166666666667 3.33333333333E2 2.38095238095E2 3.33333333333E2 7.57575757576E2
.253113553114 1.16666666667 7.09215686275 54.9711779449 529.124242424
6192.12318841 86580.2531136 1425517.16667 27298231.0678 601580873.901
15116315767.1 429614643062. 1.37116552051E13 4.88332318973E14 1.92965793419E16
8.41693047575E17 4.03380718541E19 2.11507486381E21 1.20866265223E23 7.50086674608E24 }
%%HP: T(3)A(D)F(.);
\<< { } 2. 50.
FOR n n IBERNOULLI \>NUM + 2.
STEP 'BNL' STO
\>>
This produces the list with the first twentyfive evenorder Bernoulli numbers above. The Zeta program above uses only the first four (n=8), which is suffices for the examples below. However complex arguments will require more terms as their modules grow larger, as we'll see in other examples.
Values of Zeta(s) for some real and complex arguments:
s  Zeta(s)  time(s)  exact value
+++
8.  1.00407735620  0.487  pi^8/9450
7.  1.00834927741  0.489 
6.  1.01734306199  0.493  pi^6/945
5.  1.03692775515  0.484 
4.  1.08232323370  0.488  pi^4/90
3.  1.20205690316  0.481 
2.  1.64493406685  0.479  pi^2/6
1.  9.99999999999E499  0.476  infinite
0.5  1.46035450882  0.701 
0  0.500000000000  0.388  1/2
0.5  0.207886224977  0.785 
1.  8.33333333331  0.521  1/12
2  4.00788420761E15  0.523  0
3  8.33333333323  0.535  1/120
4  2.10179397707E15  0.537  0
5  3.96825396829E3  0.543  1/252
6  2.3297903572E15  0.537  0
7  4.16666666665E3  0.540  1/240
(1.,2.)  (0.598165569758,  3.512 
 0.351854745216)  
(3.,4.)  (0.890554906967,  3.375 
 8.07594542632E3)  
(5.,6.)  (1.42843389347,  3.504 
 0.184238641590)  
Let us experiment with different numbers of zeta terms (n) and B_{n} terms (m/2) of the EulerMacLaurin summation for a few arguments in the positive complex halfplane and observe their influence on the overall accuracy using the program below.
%%HP: T(3)A(D)F(.);
\<< 0. \> s n m b
\<< '\GS(k=1.,n1.,k^s)' EVAL n 1. s  ^ s 1.  / + n s NEG ^ 2. / + BNL 1. 2. m
FOR i GETI s i + 1.  GAMMA * n 1. i  s  ^ * s GAMMA / i ! / 'b' STO+ 2.
STEP b NIP NIP +
\>>
\>>
s  n  m  Zeta(s)  time(s)
++++
2  5  4  1.64493377778  0.330
2  5  8  1.64493406547  0.438
2  7  4  1.64493403873  0.371
2  7  8  1.64493406681  0.483
2  8  10  1.64493406685  0.521
2  10  8  1.64493406685  0.464
0.5  4  4  1.46035496134  0.394
0.5  4  6  1.46035444845  0.491
0.5  4  8  1.46035451114  0.586
0.5  8  10  1.46035450881  0.743
0.5  10  8  1.46035450881  0.684
(0.5,18)  6  8  (2.32919185455,0.188602491388)  3.545
(0.5,18)  12  16  (2.32915487303,0.188866005798)  7.075
(0.5,18)  20  20  (2.32915487307,0.188866005802)  9.247
(2.5,100)  12  16  (1.02952895989,5.76320966153E2)  6.753
(2.5,100)  25  22  (1.13221432860,3.98766176236E2)  9.704
(2.5,100)  36  28  (1.13221432259,3.98766168104E2)  12.049
(0.5,14.13472514)  20  16  (2.185000000E10,1.35802000E9)  12.049
Here is the equivalent HP71B BASIC program. Because GAMMA() on the HP71B doesn't handle complex arguments, it will work only on the real domain.
10 DESTROY ALL @ OPTION BASE 1
12 REAL S,Z,B(15)
14 INTEGER I,N,M
16 INPUT S,N,M
18 FOR I=1 TO 15
20 READ B(I)
22 NEXT I
24 Z=0
26 FOR I=1 TO N1
28 Z=Z+I^(S)
30 NEXT I
32 Z=Z+N^(1S)/(S1)+N^(S)/2
34 FOR I=2 TO M STEP 2
36 Z=Z+B(I DIV 2)*GAMMA(S+I1)*N^(SI+1)/(GAMMA(S)*FACT(I))
38 NEXT I
40 PRINT Z
42 END
44 DATA 1/6,1/30,1/42,1/30,5/66,691/2730,7/6,3617/510,43867/798,174611/33
46 DATA 854513/138,236364091/2730,855303/6,23749461029/870,601580873.901
>run
? 2,5,4
1.64493377777
>run
? 0.5,4,4
1.46035496134
>run
? 2,5,8
1.64493406546
>run
? 0.5,4,8
1.46035451114
The following program HP71B program will accept both real and complex arguments:
10 DESTROY ALL @ OPTION BASE 1
12 COMPLEX P,S,Z
14 REAL B(15)
16 INTEGER I,N,M
18 INPUT S,N,M
20 FOR I=1 TO 15
22 READ B(I)
24 NEXT I
26 Z=0
28 FOR I=1 TO N1 @ Z=Z+I^(S) @ NEXT I
30 Z=Z+N^(1S)/(S1)+N^(S)/2
32 P=1/(S1)
34 FOR I=2 TO M STEP 2
36 P=P*(S+I3)*(S+I2)
38 Z=Z+B(I DIV 2)*P*N^(SI+1)/FACT(I)
40 NEXT I
42 IF IMPT(S)=0 THEN PRINT REPT(Z) ELSE PRINT Z
46 END
48 DATA 1/6,1/30,1/42,1/30,5/66,691/2730,7/6,3617/510,43867/798,174611/33
50 DATA 854513/138,236364091/2730,855303/6,23749461029/870,601580873.901
>run
? 2,5,4
1.64493377777
>run
? 0.5,4,4
1.46035496134
>run
? (0.5,18),6,8
(2.32919185455,.188602491387)
Both programs require the Math Module. I haven't tested them on the real hardware, but the running time should be fast enough.
Edited: 4 Dec 2012, 4:14 p.m.
