The Museum of HP Calculators

HP Forum Archive 21

[ Return to Index | Top of Index ]

Zeta Function [HP 50g]
Message #1 Posted by Gerson W. Barbosa on 10 Nov 2012, 5:23 p.m.

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]
Message #2 Posted by Thomas Ritschel on 11 Nov 2012, 7:09 a.m.,
in response to message #1 by Gerson W. Barbosa

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]
Message #3 Posted by Gerson W. Barbosa on 11 Nov 2012, 8:57 a.m.,
in response to message #2 by Thomas Ritschel

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]
Message #4 Posted by Eddie W. Shore on 11 Nov 2012, 9:44 a.m.,
in response to message #3 by Gerson W. Barbosa

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]
Message #5 Posted by Eddie W. Shore on 11 Nov 2012, 10:52 a.m.,
in response to message #4 by Eddie W. Shore

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]
Message #6 Posted by Gerson W. Barbosa on 11 Nov 2012, 10:53 a.m.,
in response to message #4 by Eddie W. Shore

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]
Message #7 Posted by Ángel Martin on 12 Nov 2012, 8:56 a.m.,
in response to message #6 by Gerson W. Barbosa

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]
Message #8 Posted by Thomas Ritschel on 12 Nov 2012, 2:46 p.m.,
in response to message #3 by Gerson W. Barbosa

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]
Message #9 Posted by Gerson W. Barbosa on 12 Nov 2012, 4:43 p.m.,
in response to message #8 by Thomas Ritschel

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]
Message #10 Posted by Thomas Ritschel on 14 Nov 2012, 6:59 a.m.,
in response to message #9 by Gerson W. Barbosa

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.

                                    
Re: Zeta Function [HP 50g]
Message #11 Posted by Gerson W. Barbosa on 14 Nov 2012, 9:27 p.m.,
in response to message #10 by Thomas Ritschel

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]
Message #12 Posted by Gerson W. Barbosa on 15 Nov 2012, 8:47 a.m.,
in response to message #11 by Gerson W. Barbosa

%%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]
Message #13 Posted by Gerson W. Barbosa on 19 Nov 2012, 10:11 a.m.,
in response to message #12 by Gerson W. Barbosa

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:

                  
OT (not: Zeta Function [HP 50g])
Message #14 Posted by Luiz C. Vieira (Brazil) on 14 Nov 2012, 2:32 p.m.,
in response to message #3 by Gerson W. Barbosa

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])
Message #15 Posted by Gerson W. Barbosa on 14 Nov 2012, 3:29 p.m.,
in response to message #14 by Luiz C. Vieira (Brazil)

You've got mail!

[ Return to Index | Top of Index ]

Go back to the main exhibit hall