Riemann's Zeta Function (HP 50g)



#22

This is a follow-up of a recent thread on the Zeta function . On the occasion Thomas Ritschel wrote an RPL program based on Euler-MacLaurin 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.,n-1.,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.33333333333E-2 2.38095238095E-2 -3.33333333333E-2 7.57575757576E-2
-.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 twenty-five even-order 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.00788420761E-15 | 0.523 | 0
-3 | 8.33333333323 | 0.535 | 1/120
-4 | -2.10179397707E-15 | 0.537 | 0
-5 | -3.96825396829E-3 | 0.543 | -1/252
-6 | -2.3297903572E-15 | 0.537 | 0
-7 | 4.16666666665E-3 | 0.540 | 1/240
(1.,2.) | (0.598165569758, | 3.512 |
| -0.351854745216) | |
(3.,4.) | (0.890554906967, | 3.375 |
| -8.07594542632E-3) | |
(-5.,6.) | (1.42843389347, | 3.504 |
| 0.184238641590) | |

Let us experiment with different numbers of zeta terms (n) and Bn terms (m/2) of the Euler-MacLaurin summation for a few arguments in the positive complex half-plane 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.,n-1.,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.76320966153E-2) | 6.753
(2.5,-100) | 25 | 22 | (1.13221432860,3.98766176236E-2) | 9.704
(2.5,-100) | 36 | 28 | (1.13221432259,3.98766168104E-2) | 12.049
(0.5,14.13472514) | 20 | 16 | (2.185000000E-10,-1.35802000E-9) | 12.049

Here is the equivalent HP-71B BASIC program. Because GAMMA() on the HP-71B 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 N-1
28 Z=Z+I^(-S)
30 NEXT I
32 Z=Z+N^(1-S)/(S-1)+N^(-S)/2
34 FOR I=2 TO M STEP 2
36 Z=Z+B(I DIV 2)*GAMMA(S+I-1)*N^(-S-I+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 HP-71B 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 N-1 @ Z=Z+I^(-S) @ NEXT I
30 Z=Z+N^(1-S)/(S-1)+N^(-S)/2
32 P=1/(S-1)
34 FOR I=2 TO M STEP 2
36 P=P*(S+I-3)*(S+I-2)
38 Z=Z+B(I DIV 2)*P*N^(-S-I+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.


#23

Interesting work.

Has anyone here tried the algorithms from Unified algorithms for polylogarithm, L-series, and zeta variants?


- Pauli

#24

Hi, Gerson:

Interesting work, I'll try and get some free time one of these next days (4 vacational days ahead) to try your routine, but in the meantime some well-meant, efficiency-related comments re your HP-71B versions of the code (not criticism, of course):

1) line 16 INTEGER I,N,M

INTEGER variables are slower to process and operate with than REAL variables in HP-71B BASIC, so unless memory must be saved at all costs (and with escalar, non-array variables that's never the case) using INTEGER instead of REAL is counterproductive in terms of speed.

2) lines 20 FOR I=1 TO 15, 22 READ B(I), 24 NEXT I

A FOR-NEXT loop isn't necessary to read the elements of an array in HP-71B BASIC, a simple 20 READ B will do and it will be both neater and much faster.

3) in the FOR-NEXT loops that follow, there are many invariant expressions within the loops. They should be computed and saved into variables before the loop, then use the variable inside the loop. It will be faster in general.

4) line 42 IF IMPT(S)=0 THEN PRINT REPT(Z) ELSE PRINT Z

The comparison to 0 isn't necessary, it would be more efficient and faster like this:

42 IF IMPT(S) THEN PRINT Z ELSE PRINT REPT(Z)

5) line 46 END

Not actually needed. Same with line 26 Z=0, which isn't needed after the DESTROY ALL that precedes it.

Best regards from V.



#25

Hello Valentin,

Thank you very much for your suggestions. I have only skin deep knowledge of the HP-71B system. I decided to use it before writing the HP 50g program because testing algorithms is easier in BASIC than in RPL, at least for me. I have an HP-71B with MATH ROM, 32K module and card reader, but I haven't tested the programs on it yet. I suspect the two calls to the Gamma function in the last loop takes longer than the alternative two multiplication operations, even for a large number of iteration, but I choose it for compactness. Also, pre-computed Bernoulli numbers is mandatory on the HP-71B, but when only four of them are needed perhaps the built-in BERNOULLI function on the HP 50g would be faster. Anyway, half a second running time for arguments in the limited range most users would be interested in fine for me. That's about the time trigs took on the HP-35 :-)
I have the book in the reference since 2009, but I never read past the first chapter (The math is too complicated for me).

Any improvement you want to try and share, such as automatic choice of parameter for the maximum possible range, is much appreciated.

Best regards,

Gerson.


#26

Hello Gerson,

thanks again for sharing your work. It's interesting to see that you got the algorithm much faster by precomputing the Bernoulli numbers.

Meanwhile I tried the Borwein algorithm. There is also a related paper by Xavier Gourdon and Pascal Sebah.


For some initial tests I prepared this little Python script:

#!/usr/bin/python
# computes the zeta function for complex arguments
# using Borwein's "Algorithm 2"

import math, cmath, sys

s = complex(float(sys.argv[1]),float(sys.argv[2]))

def d_k(k,n):
y = 0.0
for i in range(0,k+1):
y += math.factorial(n+i-1)*(4.0**i)/(math.factorial(n-i)*math.factorial(2*i))
y *= n
return y

def zeta2(s,n):
y = 0.0
for k in range(0,n):
y += ((-1)**k) * (d_k(k,n)-d_k(n,n)) / (k+1)**s
y *= (-1)/(d_k(n,n) * (1-2**(1-s)))
return y

for n in range(1,25):
print n, zeta2(s,n)

It turns out that for 12 digits precision one needs about n=12-15.
This is in accordance with n=1.3d+0.9t suggested by Gourdon and Sebah ("d" is the number of digits, "t" is the absolute value of the imaginary part of the argument).

On the HP50g the same algorithm as given above takes roughly half a minute for n=14, which is quite slow. Perhaps, one could try some optimization by unrolling the double loop, or at least reducing the number of "d_k" evaluations.

Borwein presents in his paper also a "simpler" algorithm (dubbed "algorithm 3") using binomial coefficents. But there are still double sums, which makes it quite slow for n>10.

Perhaps I should consider this as a challenge to dive a little deeper into SysRPL, in which I have almost no experience so far...

Best regards,

Thomas.


#27

The 34S uses the Borwein algorithm. We found that the error bounds he provides are too tight and more iterations are actually required. However it is still the fastest algorithm I uncovered at the time. The other problem is this algorithm takes O(y) iterations to converge for a complex argument x + iy. This is why the 34S doesn't support complex zeta even though I did code it.


- Pauli


#28

Quote:
This is why the 34S doesn't support complex zeta even though I did code it.

This is the answer to the question I was going to ask you. Thank you!

Gerson.


#29

Knowing that the function will execute in reasonable time is kind of nice. I couldn't guarantee this with complex zeta unfortunately.

- Pauli


#30

Yes but even if slow it's a pleasure to obtain those results. The 41Z has it implemented in ZZETA, it works like a charm on the CL.


#31

Someone will try to evaluate 3 + i 1e25...


- Pauli


#32

easy: "DATA ERROR"
obviously ZZETA is not up to this challenge :-)

http://www.wolframalpha.com/input/?i=polylog%283+%2B+i+1e25%29


#33

Neither is Zeta above, which cannnot properly handle 0.5+18i :-)


#34

ZZETA obtains in a blitz (on the CL) for that argument:

18, ENTER, .5, ZZETA --> 2,329-J0,189

and using FIX9 settings:

Im = -0,188866015

Re = 2,329154846

not bad for a venerable platform.


#35

The first HP 50g program above, with fix parameters n=10 and m=8, gets the first 5 or 6 digits right:

(2.32915540368,-0.188867886498)
This takes about 4 seconds. As Im(s) grows larger, the accuracy degrades significantly.
The actual 12-digit result is
(2.32915487305,-0.18886600580),
which means your best HP-41 result is really not bad:
(2.329154846,-0.188866015)


#36

Using n=30 the Borwein algorithm (see the program in message #18) yields:

zeta(0.5+18i) -> (2.32915487306,-0.188866005808)

This took about 4 seconds on a HP 49g+.


#37

That's pretty good! I think it will be the only Zeta program in my 50g. Thanks again!

#38

Hello Thomas,

Quote:
thanks again for sharing your work.

It's me who should thank you again. This is essentially the work you have started. I only recoded it following the text and examples in my reference.

Quote:
It's interesting to see that you got the algorithm much faster by precomputing the Bernoulli numbers.

Probably the idea wouldn't have occurred to me if I hadn't decided to test the algorithm on the HP-71B first. When doing it I noticed the computation of B2, B4,...B10 alone would take about one minute.

Quote:
Meanwhile I tried the Borwein algorithm. There is also a related paper by Xavier Gourdon and Pascal Sebah.

That's the way to go. This is the same method used in SandMath, I'm told. Please share your HP 50g when it's finished.


Quote:
Perhaps I should consider this as a challenge to dive a little deeper into SysRPL, in which I have almost no experience so far...

Besides speed another advantage of SysRPL is the ability of doing all calculations with 15 significant digits, I think. I don't know anything about SysRPL, though.

Best regards,

Gerson.


#39

Quote:
Besides speed another advantage of SysRPL is the ability of doing all calculations with 15 significant digits, I think.

SYSRPL has the long real and long complex types which do carry extra digits -- I think fifteen but could be mistaken.


- Pauli

#40

Hello Gerson,

here comes a (UserRPL) implementation of the Borwein algorithm for the HP 50g. Actually it's the Gourdon-Sebah variant ("Proposition 1" in their paper). Note that their definition of "d_k" is slightly different to Borwein's.

DK:

%%HP: T(3)A(R)F(.);
\<< \-> k n
\<< 0 k n
FOR i n i + 1 - !
4 i ^ * n i - ! 2 i *
! * / +
NEXT n *
\>>
\>>


Zeta:

%%HP: T(3)A(R)F(.);
\<< \-> s n
\<< 0 1 n
FOR k -1 k 1 - ^
k n DK * k s ^ / +
NEXT 1 0 n DK 1 2
1 s - ^ - * / *
\>> EVAL \->NUM
\>>

The above version is more or less just a proof of concept.
For n=10 or more iterations it is rather slow.

However, it can be improved by precomputing the "d_k".
The following routine generates the list of "d_k" for a given n. Just enter the number of iterations, e.g. n=20, and run "DKI":

DKI:

%%HP: T(3)A(R)F(.);
\<< \-> n
\<< { } 0 n
FOR i i n DK EVAL
+
NEXT
\>> 'DKL' STO
\>>

The list is stored in "DKL" (for "d_k list") and is then used by the following program. It automatically takes the number of iterations from the size of the "DKL" list.

Zeta2:

%%HP: T(3)A(R)F(.);
\<< PUSH -105 SF DKL
SIZE 1 - \-> s n
\<< DKL 1 GETI 1 2 1
s - ^ - * INV NEG 3
ROLLD 1 n
FOR k GETI -1 k ^
* k s ^ / 3 ROLLD
NEXT DROP2 2 n
FOR i +
NEXT *
\>> EVAL \->NUM POP
\>>

It should work for real and complex arguments and is reasonable fast. I did no timings, but it "feels" like less than 2 seconds for 20 iterations. Most probably the EVAL and \->NUM steps could be removed.

Best regards,

Thomas

#41

Hello Valentin,

Quoting myself above:

Quote:
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)

Dover Publications, Inc. has kindly granted me permission to reproduce the following excerpt for use by forum members of this museum:

http://i918.photobucket.com/albums/ad30/gwbarbosa/Temp/Zeta_E-McL_Sum_1.png

http://i918.photobucket.com/albums/ad30/gwbarbosa/Temp/Zeta_E-McL_Sum_2.png

http://i918.photobucket.com/albums/ad30/gwbarbosa/Temp/Zeta_E-McL_Sum_3.png

http://i918.photobucket.com/albums/ad30/gwbarbosa/Temp/Zeta_E-McL_Sum_4.png

http://i918.photobucket.com/albums/ad30/gwbarbosa/Temp/Zeta_E-McL_Sum_5.png

Riemann's Zeta Function, by H. M. Edwards. Dover Publications (2001), ISBN 10: 0486417409, pages 114 through 118.

I think this book might be useful to you and other readers who have stronger math backgrounds than I do.

Best regards,

Gerson.


#42

Implementation of the above method on the HP-28S and HP-48 calculators:

Zeta (HP-48)
%%HP: T(3)A(D)F(.);
\<< DUP 1 - INV 10 8 0 \-> s p n m b
\<< 0 1 n 1 -
FOR k k s NEG ^ +
NEXT n 1 s - ^ s 1 - / + n s NEG ^ 2 / + BNL 1 2 m
FOR i s i + 3 - DUP 1 + * 'p' STO* GETI p * n 1 i - s - ^ * i ! / 'b' STO+ 2
STEP b ROT ROT DROP2 +
\>>
\>>

Zeta (HP-28S)
\<< DUP 1 - INV 10 8 0 \-> s p n m b
\<< 0 1 n 1 -
FOR k k s NEG ^ +
NEXT n 1 s - ^ s 1 - / + n s NEG ^ 2 / + BNL 1 2 m
FOR i s i + 3 - DUP 1 + * p * 'p' STO GETI p * n 1 i - s - ^ * i FACT / b + 'b' STO 2
STEP b ROT ROT DROP2 +
\>>
\>>

BNL
{ .166666666667 -3.33333333333E-2 2.38095238095E-2 -3.33333333333E-2 7.57575757576E-2
-.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 }

Re(s) >= -1; small Im(s)

For larger complex arguments, edit the program and change the n and m parameters as needed (up to m=50). The default values (n=10, m=8) should give maximum accuracy for real arguments (x>=-1).

----------------------------------------------------

P.S.: Running time on the HP-28S

(0.5 18) Zeta  --> (2.32915540368, -0.188867886498)  in 4.8610 seconds  (average of three timings using # 4554d SYSEVAL)
By replacing
FOR i s i + 3 - DUP 1 + * p * 'p' STO
with
FOR i s i + 3 - DUP SQ + p * 'p' STO
the running time gets down to 4.7005 seconds.

Other such minor improvements can be done, but what would probably halve the running time would be using the stack only.


Edited: 19 Dec 2012, 12:17 p.m.


Possibly Related Threads…
Thread Author Replies Views Last Post
  HP50g: Writing a function that returns a function Chris de Castro 2 2,162 12-10-2013, 06:49 PM
Last Post: Han
  IFERR function on HP Prime Mic 2 1,863 12-02-2013, 01:33 AM
Last Post: cyrille de Brébisson
  HP Prime: Dirichlet's eta function recognized but not numerically evaluated Helge Gabert 0 1,169 11-16-2013, 03:41 PM
Last Post: Helge Gabert
  Possible bug with sqrt function in the HP prime Michael de Estrada 6 2,568 11-15-2013, 12:49 PM
Last Post: Michael de Estrada
  HP-41 MCODE: The Last Function - at last! Ángel Martin 0 1,101 11-08-2013, 05:11 AM
Last Post: Ángel Martin
  HP Prime 'where' function bluesun08 11 3,485 10-29-2013, 06:56 PM
Last Post: Joe Horn
  HP Prime - Defining a function bluesun08 5 2,347 10-23-2013, 02:43 PM
Last Post: Han
  HP PRIME Function parameters steindid 2 1,473 10-11-2013, 10:20 AM
Last Post: steindid
  HP Prime can't find how to use partfrac function ! dg1969 3 1,679 10-04-2013, 09:25 PM
Last Post: Joe Horn
  HP Prime function APP - Strange limitation ! :o( dg1969 2 1,629 10-04-2013, 12:10 PM
Last Post: dg1969

Forum Jump: