HP Forums

Full Version: Riemann's Zeta Function (HP-28S)
You're currently viewing a stripped down version of our content. View the full version with proper formatting.

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.

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

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.

On my 28S now :-)


- Pauli

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.

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

- Pauli

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.

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.

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.