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: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.%%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

\>>

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 B_{n} 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. *