Riemann's Zeta Function (HP-28S)

[Gerson W. Barbosa]

For the ones who still cherish the good old HP-28S after all those years :-)

Not at all optimized, but it runs reasonably fast. Please test for yourself the range where it is accurate enough for your needs.

[Eduardo Duenez]

Thanks. Big Fan both of the Zeta function and the 28S, which was the first HP calculator I ever owned (then sold for a pittance, then wonderfully re-acquired perchance from Julian of Spain via this Forum).

Eduardo

[Gerson W. Barbosa]

The HP-28S was my second HP calculator, replacing a HP-15C (now I have both back, the HP-15C being my original one).
I may post the text listing later (The only change is the RAD instruction being placed in the beginning of the last IF-THEN clause).
The major drawback of the HP-28S is the lack of a PC link. I occasionally lost everything in my original 28S due to a bad battery contacly, fortunately my new unit don't have this problem. It was an advanced scientific calculator then and still is today: it can even plot the graph in the end of the archived thread below, albeit its smaller screen :-)

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv021.cgi?read=233799#233799

There is an older HP-28S version here:

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv021.cgi?read=235199

Best regards,

Gerson.

[Paul Dale]

On my 28S now :-)

- Pauli

[Gerson W. Barbosa]

The Gamma function program gives 6 or 7 digits of accuracy for x = 2 and 11 for x = 4 and greater (positive arguments only). I used the approximation involving the sinh(x) function in Viktor Toth's page. Additional correction terms were obtained by comparing the sinh(x) expansion with this series:

series[e^(2*(1/(12*x^2)-1/(360x^4)+1/(1260*x^6)-1/(1680*x^8)+1/(1188*x^10)-691/(360360*x^12)+1/(156*x^14)-3617/(122400*x^16)+43867/(244188*x^18)-174611/(125400*x^20)+77683/(5796*x^22)))]

W|A series expansion at x = oo

Gerson.

[Paul Dale]

A question to the RPL users: why do people use lower case for local variables?

- Pauli

[Gerson W. Barbosa]

Quoting from the HP-28S Owner's Manual, page 80:

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

It's useful to follow some convention to distinguish your local variables from your ordinary or "global" variables. This manual uses lower-case letters to distinguish local variables.

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

Gerson.

[Gerson W. Barbosa]

Here are the listings, in case the picture becomes unavailable:

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

HP-28S 

Zeta:

\<< RCLF 36 CF SWAP 1

CF DUP RE -1 <

  IF

  THEN 1 SF NEG 1 +

  END DUP 1 - INV

OVER NEG 10 8 0 \-> s

p q n m b

  \<< 0 1 n 1 -

    FOR k k q ^ +

    NEXT n 1 s - ^ s 

1 - / + n q ^ 2 / + 

BNL 1 2 m

    FOR i s i + 3 -

DUP SQ +  p * 'p' STO 

GETI p * n 1 i - s - 

^ * i FACT / b + 'b'

STO 2

    STEP b ROT ROT

DROP2 + 1 FS?

    IF

    THEN 1 CF RAD 1

s - DUP 2 / \pi * SIN

2 ROT ^ * * s DUP

Gamma \pi ROT NEG ^ *

*

    END

  \>> SWAP STOF

\>>

Gamma:

\<< DUP SQ \-> z z2

  \<< z DUP DUP '(z2*(

(z2*((54486432000*z2

-42291849600)*z2+

66097231200)-

157264451760)*z2+

538162408630)-

2.51535889068E12)/(

4.413400992E13*SQ(SQ

(SQ(z2))))' EVAL z 

DUP INV SINH * + \v/ * 

e / SWAP ^ \pi DUP + 

ROT / \v/ * \->NUM

  \>>

\>>

BNL:

{ .166666666667                 ;  1/6                       

-3.33333333333E-2               ; -1/30                            

2.38095238095E-2                ;  1/42

-3.33333333333E-2               ; -1/30

7.57575757576E-2                ;  5/66

-.253113553114                  ; -691/2730

1.16666666667                   ;  7/6 

-7.09215686275                  ; -3617/510 

54.9711779449                   ;  43867/798

-529.124242424 }                ; -174611/330 

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

HP-48G/GX

Zeta:

%%HP: T(3)A(D)F(,);

\<< RCLF SWAP 1 CF

DUP RE -1 <

  IF

  THEN 1 SF NEG 1 +

  END DUP 1 - INV

OVER NEG 10 8 0 \-> s

p q n m b

  \<< 0 1 n 1 -

    FOR k k q ^ +

    NEXT n 1 s - ^

s 1 - / + n q ^ 2 /

+ BNL 1 2 m

    FOR i s i + 3 -

DUP SQ + 'p' STO*

GETI p * n 1 i - s

- ^ * i ! / 'b'

STO+ 2

    STEP b ROT ROT

DROP2 + 1 FS?

    IF

    THEN 1 CF RAD 1

s - DUP 2 / \pi * SIN

2 ROT ^ * * s DUP

Gamma \pi ROT NEG ^ *

*

    END

  \>> SWAP STOF

\>>

Gamma:

%%HP: T(3)A(D)F(,);

\<< DUP SQ \-> z z2

  \<< z DUP DUP '(z2*

((z2*((54486432000*

z2-42291849600)*z2+

66097231200)-

157264451760)*z2+

538162408630)-

2,51535889068E12)/(

4,413400992E13*SQ(

SQ(SQ(z2))))' EVAL

z DUP INV SINH * +

\v/ * e / SWAP ^ \pi

DUP + ROT / \v/ *

  \>>

\>>

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 }

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

Edited: 29 Jan 2013, 9:39 p.m.

[Gerson W. Barbosa]

Because of finite precision -- sine(x*pi/2) never evaluates to zero for even x as expected -- the program doesn't return the trivial zeroes properly. For instance,

 Zeta( -2.000)  -->   4.00788420477E-15

 Zeta( -4.000)  -->  -2.10179397700E-15

 Zeta( -6.000)  -->   2.32972903567E-15

 ...

 Zeta(-30.000)  -->  -4.49028406857E-3

 ...

 Zeta(-59.999)  -->   5.33788222631E30

 Zeta(-60.000)  -->  -2.11263214013E22  (way off!)

 Zeta(-60.001)  -->  -5.36211534566E30

In order to prevent this from happening, insert this optional patch between SIN and 2, in the main program:

...

s - DUP 2 / \pi * SIN

PATCH 2 ROT ^ * * s 

...

PATCH:

%%HP: T(3)A(D)F(,);

\<< OVER DUP IM NOT

  \<< RE DUP 2 MOD

NOT NOT * SIGN NEG

*

  \>>

  \<< DROP

  \>> IFTE

\>>

Example:

 Zeta(-59.999)  -->   5.33788222631E30

 Zeta(-60.000)  -->   0.00000000000

 Zeta(-60.001)  -->  -5.36211534566E30

Also, in order to avoid wrong answers due to underflow and overflow, make sure flags 57 and 58 (-20 and -21 on the HP-48) are set.