Zeta Function [HP 50g]



#16

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.


#17

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, ...}


#18

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.


#19

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


#20

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

#21

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.


#22

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.

#23

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


#24

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.


#25

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.


#26

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.


#27

%%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


#28

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:

#29

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)


#30

You've got mail!


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

Forum Jump: