Riemann's Zeta Function (HP-28S) « Next Oldest | Next Newest »

 ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 01-25-2013, 09:28 PM 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 Member Posts: 56 Threads: 10 Joined: May 2006 01-28-2013, 06:48 PM 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 Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 01-28-2013, 07:25 PM 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 :-) There is an older HP-28S version here: Best regards, Gerson. Paul Dale Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 01-28-2013, 07:38 PM On my 28S now :-) - Pauli ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 01-28-2013, 08:14 PM 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 Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 01-29-2013, 02:21 AM A question to the RPL users: why do people use lower case for local variables? - Pauli ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 01-29-2013, 08:37 AM 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 Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 01-29-2013, 09:28 PM 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 Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 02-03-2013, 03:23 PM 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 1,037 12-10-2013, 06:49 PM Last Post: Han IFERR function on HP Prime Mic 2 927 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 598 11-16-2013, 03:41 PM Last Post: Helge Gabert Possible bug with sqrt function in the HP prime Michael de Estrada 6 1,182 11-15-2013, 12:49 PM Last Post: Michael de Estrada HP-41 MCODE: The Last Function - at last! Ángel Martin 0 549 11-08-2013, 05:11 AM Last Post: Ángel Martin HP Prime 'where' function bluesun08 11 1,833 10-29-2013, 06:56 PM Last Post: Joe Horn HP Prime - Defining a function bluesun08 5 1,125 10-23-2013, 02:43 PM Last Post: Han HP PRIME Function parameters steindid 2 745 10-11-2013, 10:20 AM Last Post: steindid HP Prime can't find how to use partfrac function ! dg1969 3 798 10-04-2013, 09:25 PM Last Post: Joe Horn HP Prime function APP - Strange limitation ! :o( dg1969 2 820 10-04-2013, 12:10 PM Last Post: dg1969

Forum Jump: