▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
▼
Posts: 17
Threads: 1
Joined: Nov 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, ...}
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
▼
Posts: 764
Threads: 118
Joined: Aug 2007
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
▼
Posts: 764
Threads: 118
Joined: Aug 2007
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
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
▼
Posts: 1,253
Threads: 117
Joined: Nov 2005
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.
Posts: 17
Threads: 1
Joined: Nov 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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
▼
Posts: 17
Threads: 1
Joined: Nov 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.
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
%%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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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:
Posts: 591
Threads: 16
Joined: Feb 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)
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
|