 Zeta Function [HP 50g] - Printable Version +- HP Forums (https://archived.hpcalc.org/museumforum) +-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html) +--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html) +--- Thread: Zeta Function [HP 50g] (/thread-233799.html) Zeta Function [HP 50g] - Gerson W. Barbosa - 11-10-2012 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. Re: Zeta Function [HP 50g] - Thomas Ritschel - 11-11-2012 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, ...}  Re: Zeta Function [HP 50g] - Gerson W. Barbosa - 11-11-2012 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. Re: Zeta Function [HP 50g] - Eddie W. Shore - 11-11-2012 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 Re: Zeta Function [HP 50g] - Eddie W. Shore - 11-11-2012 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  Re: Zeta Function [HP 50g] - Gerson W. Barbosa - 11-11-2012 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. Re: Zeta Function [HP 50g] - Ángel Martin - 11-12-2012 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. Re: Zeta Function [HP 50g] - Thomas Ritschel - 11-12-2012 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 Re: Zeta Function [HP 50g] - Gerson W. Barbosa - 11-12-2012 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. Re: Zeta Function [HP 50g] - Thomas Ritschel - 11-14-2012 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. OT (not: Zeta Function [HP 50g]) - Luiz C. Vieira (Brazil) - 11-14-2012 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) Re: OT (not: Zeta Function [HP 50g]) - Gerson W. Barbosa - 11-14-2012 You've got mail! Re: Zeta Function [HP 50g] - Gerson W. Barbosa - 11-14-2012 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.  Re: Zeta Function [HP 50g] - Gerson W. Barbosa - 11-15-2012 %%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  Re: Zeta Function [HP 50g] - Gerson W. Barbosa - 11-19-2012 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: 