Riemann's Zeta Function (HP-28S)



#6

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.


#7

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


#8

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.

#9

On my 28S now :-)


- Pauli


#10

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.

#11

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

- Pauli


#12

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.

#13

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.


#14

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.


Possibly Related Threads…
Thread Author Replies Views Last Post
  HP50g: Writing a function that returns a function Chris de Castro 2 2,166 12-10-2013, 06:49 PM
Last Post: Han
  IFERR function on HP Prime Mic 2 1,866 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,569 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: