Zeta Function [HP 50g] Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 11-10-2012, 05:23 PM Here are two HP 50g programs to compute Zeta(x). The first one will accept real or complex arguments (no extended domain, however), but it's too slow for small arguments (Re(x) < 4). The second one is faster but will accept only a limited number of integers. %%HP: T(3)A(D)F(.); \<< 11.4 OVER RE LN / 3. - EXP IP 1. SWAP 1. 1. ROT START NEXTPRIME PICK3 OVER SWAP ^ INV NEG 1. + ROT * SWAP NEXT DROP INV NIP \>> x | Zeta(x) | time(s) ---------+--------------------+--------- 3 | 1.2020569029 | 89.17 3.5 | 1.12673386727 | 21.75 4 | 1.08232323369 | 8.49 5 | 1.03692775514 | 2.57 6 | 1.01734306198 | 1.20 7 | 1.00834927738 | 0.74 (3,5) | (0.912526588591, | 178.42 | 5.08428709213E-2) | (8,9) | (1.00377981783, | 1.04 | 2.42567436216E-4) | Re(x) > 1 (but it will take infinite time when x approaches 1) %%HP: T(3)A(D)F(.); \<< 1. - DUP 1. SWAP PSI SWAP ! / DUP SIGN * \>> x | Zeta(x) | time(s) ---------+--------------------+--------- 2 | 1.64493406685 | 0.27 3 | 1.20205690316 | 0.27 4 | 1.08232323371 | 0.27 5 | 1.03692775514 | 0.28 6 | 1.01734306198 | 0.28 7 | 1.00834927738 | 0.28 x = {2, 3, 4,... 11}  Edited: 10 Nov 2012, 5:31 p.m. Thomas Ritschel Junior Member Posts: 17 Threads: 1 Joined: Nov 2012 11-11-2012, 07:09 AM Thanks for sharing your programs for computing the zeta function, Gerson! One could slightly modify the second one so that it would also work for larger integer values (but starting from x=12 it's limited to even arguments): %%HP: T(3)A(D)F(.); \<< 1 - DUP 1 SWAP PSI ->NUM SWAP ! / DUP SIGN * \>> x | Zeta(x) ---------+-------------------- 8 | 1.0040773562 9 | 1.00200839282 10 | 1.00099457513 11 | 1.00049418861 12 | 1.00024608655 14 | 1.00006124814 16 | 1.00001528226 18 | 1.00000381729 20 | 1.00000095396 x = {2, 3, 4, ... 11, 12, 14, 16, ...}  Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 11-11-2012, 08:57 AM Hallo Thomas, I remember my second program used to work for larger even numbers, but when I edited in approximate mode it became limited to 12. I had changed it because it would always force exact mode on. The ->NUM instruction did the trick. Thanks! The first program is just an implementation of the Euler's Product Form, as you can see. I hope someone comes up with a really fast implementation of Riemann's Zeta Function for the HP-50g. This would be quite welcome! Regards, Gerson. Eddie W. Shore Posting Freak Posts: 764 Threads: 118 Joined: Aug 2007 11-11-2012, 09:44 AM Gerson, I was also able to use larger integers with the second program - the exact mode does the trick. This is awesome. Here is what I have - trying to take the direct route (but I capped the upper limit to reduce calculation time). The approximation gets more accurate as X gets higher (7 decimal places at x=4, almost full machine accuracy at x=12), but is pathetic when x is small (x < = 2). << PUSH -105 SF - > Z 'Sigma(X=1,500,1/X^Z)' EVAL POP>>  Eddie Eddie W. Shore Posting Freak Posts: 764 Threads: 118 Joined: Aug 2007 11-11-2012, 10:52 AM Here is version where: 1. If the argument is an integer (fractional part is 0), then Gerson Barbosa's second program is executed and Exact Mode is turned on (Flag -105 is cleared). 2. Otherwise, the sum is calculated until either display accuracy (12 decimal places) is reached or 1,500 iterations happened. This still is not so desirable for x < 2. The limits can be adjusted if you want a shorter calculation time.  << -> Z << IF Z FP 0 == THEN -105 CF Z 1 - DUP 1 SWAP PSI SWAP ! / DUP SIGN * ELSE 0 'K' STO 0 DO 'K' INCR Z NEG ^ + EVAL UNTIL K Z NEG ^ 1E-13 < 1500 K == OR END 'K' PURGE END >> >> 214.5 bytes  Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 11-11-2012, 10:53 AM Hello Eddie, The Euler Product Form requires less terms but won't perform any better for small arguments. Also those many multiplications will significantly worsen the accuracy in that range. Ángel has a complete implementation of Riemann's Zeta function on his SandMath package. I don't know how fast it is on the real HP-41, but it is fast enough on the emulator. I guess an RPL conversion of the algorithm therein would be very fast on the HP-48/49/50g. Regards, Gerson. Ángel Martin Posting Freak Posts: 1,253 Threads: 117 Joined: Nov 2005 11-12-2012, 08:56 AM There are in fact two implementations in the SandMath, one (ZETA) using a direct method in MCODE and another (ZETAX) using the Borwein algorithm with a FOCAL program. ZETAX is more accurate and faster for smaller argument; and also capable of calculating ZETA for very small arguments. All I did was adapt JM Baillard's program to the SandMath environment, so it's not my credit. The MCODE ZETA is more of a demo than a practical use. Run it with FIX 9 set and you'll see the digits flick... Edited: 12 Nov 2012, 8:59 a.m. Thomas Ritschel Junior Member Posts: 17 Threads: 1 Joined: Nov 2012 11-12-2012, 02:46 PM Hello Gerson, the german wikipedia article on the zeta function provides an Euler-McLaurin sum formula, which can easiliy programmed into the HP-50g: \<< \-> s N m \<< 0 1 N 1 - FOR n n s NEG ^ + NEXT N 1 s - ^ s 1 - / + 1 2 / N s NEG ^ * + 1 m FOR m 2 m * DUP IBERNOULLI SWAP ! / 0 2 m * 2 - FOR n s n + * NEXT N s NEG 2 m * - 1 + ^ * + NEXT ->NUM \>> \>>  Note that the above code requires 3 arguments on the stack: besides the argument of the zeta function (called "s" in the code above) also the two integers "N" and "m" which more or less control the number of terms in the sum and therefore the accuracy. For example, for s=1/2, N=4, and m=2, one will get: zeta(1/2) = -1.46035496134 Using N=6 and m=4 leads to: zeta(1/2) = -1.46035450885 I did no precise timings so far, but the two examples took no longer than about 1 or 2 seconds, respectively. The above code could be further optimized, e.g. by using the fact that the successive "Bernoulli" terms share common factors. Kind regards, Thomas Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 11-12-2012, 04:43 PM Thank you very much, Thomas! That's what I was looking for. N=8 and m=4 appears to give maximum accuracy most of the times. This is also very fast:  1.0000001 --> 10000000.5772 (1.75s) 2 --> 1.64493406684 (1.77s) (3,5) --> (.91256588997, 5.08428710748E-2) (3.42s) (.5,14.134725) --> (1.76714676673E-8, (N=14, m=7) -1.11022625159E-7) (6.43s)  Best regards, Gerson. Thomas Ritschel Junior Member Posts: 17 Threads: 1 Joined: Nov 2012 11-14-2012, 06:59 AM Hello Gerson, I replaced the nested for-loops by a recursive computation of the "Bernoulli" terms and got a nice little speed improvement (up to about 25% depending on the number of terms): %%HP: T(3)A(D)F(.); \<< \-> s N m \<< 0 1 N 1 - FOR n n s NEG ^ + NEXT N 1 s - ^ s 1 - / + 0 m 2 FOR n 2 n * IBERNOULLI + 2 n * DUP 1 - DUP 1 - DUP 1 - s + SWAP s + * SWAP / SWAP / N DUP * / * -1 STEP 2 IBERNOULLI + s N / * 1 + 2 / N s ^ / + EVAL \->NUM \>>  Note the "EVAL" at the end of the code. I noticed deviations of +/-1 in the last decimal place if the "evaluation step" is omitted. Here are some timings:  x | N | m | Zeta(x) | time(s) (old/new) -----+-----+-----+---------------+----------------- 2 | 5 | 2 | 1.64493377778 | 0.29 / 0.23 2 | 5 | 4 | 1.64493406547 | 0.67 / 0.49 2 | 5 | 6 | 1.64493406682 | 1.38 / 0.99 2 | 5 | 8 | 1.64493406684 | 2.58 / 1.88 2 | 5 | 10 | 1.64493406685 | 4.33 / 3.24  Kind regards, Thomas. Luiz C. Vieira (Brazil) Senior Member Posts: 591 Threads: 16 Joined: Feb 2012 11-14-2012, 02:32 PM Gerson, you OK? I sent you a message through the MoHPC message system. Did you get it? I'd just like to know that. Cheers. Luiz (Brazil) Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 11-14-2012, 03:29 PM You've got mail! Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 11-14-2012, 09:27 PM Hello Thomas, I've made a small modification to improve the accuracy for negative arguments. The program keeps asking for approximate mode change when the argument is non-integer. Perhaps a convergence test should be included (this is better than fixed N and m parameters). Also, PUSH and POP, as Eddie used in one of his programs, should be used to save and restore the user's flags. Best regards, Gerson. %%HP: T(3)A(R)F(,); \<< 1 CF DUP RE 0 < IF THEN 1 SF NEG 1 + END 5 8 \-> s N m \<< 0 1 N 1 - FOR n n s NEG ^ + NEXT N 1 s - ^ s 1 - / + 0 m 2 FOR n 2 n * IBERNOULLI + 2 n * DUP 1 - DUP 1 - DUP 1 - s + SWAP s + * SWAP / SWAP / N DUP * / * -1 STEP 2 IBERNOULLI + s N / * 1 + 2 / N s ^ / + 1 FS? IF THEN 1 s - DUP 2 / \pi * SIN 2 ROT ^ * * s DUP GAMMA \pi ROT NEG ^ * * END EVAL \->NUM \>> \>> -7 --> 4.16666666666E-3 1/x --> 240. -3 --> 8.33333333331E-3 1/x --> 120. -9 --> -7.57575757575E-3 1/x --> -132.  Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 11-15-2012, 08:47 AM %%HP: T(3)A(D)F(.); \<< PUSH RAD 1. CF DUP RE 0. < IF THEN 1. SF NEG 1. + END 5. 9. \-> s N m \<< 0. 1. N 1. - FOR n n s NEG ^ + NEXT N 1. s - ^ s 1. - / + 0. m 2. FOR n 2. n * IBERNOULLI -105. SF + 2. n * DUP 1. - DUP 1. - DUP 1. - s + SWAP s + * SWAP / SWAP / N DUP * / * -1. STEP 2. IBERNOULLI -105. SF + s N / * 1. + 2. / N s ^ / + 1. FS? IF THEN 1. s - DUP 2. / \pi * SIN 2. ROT ^ * * s DUP GAMMA \pi ROT NEG ^ * * END \>> POP \>> 2. --> 1.64493406685 (6.10s) -1. --> -8.33333333331E-3 (6.14s) 1/x --> -12. -3. --> 8.33333333335E-3 (6.17s) 1/x --> 120. -5. --> -3.96825396829E-3 (6.17s) 1/x --> -251.999999998 -7. --> 4.16666666665E-3 (6.16s) 1/x --> 240.000000001 -9. --> -7.5757575758E-3 (6.16s) 1/x --> -131.999999999  Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 11-19-2012, 10:11 AM The first three Riemann Zeta zeroes, (0.5 14.13), (0.5 21.02) and (0.5 25.01), can be nicely plotted on the HP 50g using Thomas Ritschel's program above and the following plot settings: « Next Oldest | Next Newest »

 Possibly Related Threads... Thread Author Replies Views Last Post HP50g: Writing a function that returns a function Chris de Castro 2 1,610 12-10-2013, 06:49 PM Last Post: Han IFERR function on HP Prime Mic 2 1,434 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 879 11-16-2013, 03:41 PM Last Post: Helge Gabert Possible bug with sqrt function in the HP prime Michael de Estrada 6 1,840 11-15-2013, 12:49 PM Last Post: Michael de Estrada HP-41 MCODE: The Last Function - at last! Ángel Martin 0 795 11-08-2013, 05:11 AM Last Post: Ángel Martin HP Prime 'where' function bluesun08 11 2,703 10-29-2013, 06:56 PM Last Post: Joe Horn HP Prime - Defining a function bluesun08 5 1,786 10-23-2013, 02:43 PM Last Post: Han HP PRIME Function parameters steindid 2 1,141 10-11-2013, 10:20 AM Last Post: steindid HP Prime can't find how to use partfrac function ! dg1969 3 1,158 10-04-2013, 09:25 PM Last Post: Joe Horn HP Prime function APP - Strange limitation ! :o( dg1969 2 1,183 10-04-2013, 12:10 PM Last Post: dg1969

Forum Jump: