Multidimensional Zeta-Functions for the HP-41

Overview
 

 1°)  Barnes Zeta Function

  a)  Series
  b)  Integral Representation

 2°)  Epstein Zeta Function
 3°)  Multiple Zeta Function

  a)  Program#1
  b)  Program#2
 
 

1°)  Barnes Zeta Function
 

     a) Series
 

-Barnes' Zeta function is defined by

    z_B( s,w | a₁ , ... , a_n ) = SUM_{k1,....,kn=0,1,2,....}    1 / ( w + a₁ k₁ + a₂ k₂ + ............... + a_n k_n ) ^s

    and    s > n   ,  w >= 0
 

-If w = 0 , the term corresponding to   k₁ = k₂ = .......... = k_n = 0  is omitted
 

Data Registers:           •  R00 = n                                             ( Registers R00 thru Rnn are to be initialized before executing "BZF" )

                                      •  R01 = a₁     •  R02 = a₂    ......................       •  Rnn = a_n
                                         Rn+1 = k₁      Rn+2 = k₂  ......................          R2n = k_n
Flags: /
Subroutines: /
 
 

 
     ( 134 bytes / SIZE 2n+1 )
 
 

   Provided  R00 = n  ,   R01 = a₁ ,  ............ , Rnn = a_n

Examples:

    •  z_B( 7 , 0.8 | 4 , 5 , 6 ) = ?

    n = 3   STO 00    4  STO 01   5  STO 02   6  STO 03

     7    ENTER^
   0.8   XEQ "BZF"  >>>>  z_B( 7 , 0.8 | 4 , 5 , 6 ) = 4.768395227           ---Execution time = 54s---

-Likewise,

    •  z_B( 7 , 1.2 | 4 , 5 , 6 ) = 0.279095669                          ---Execution time = 97s---
    •  z_B( 7 , 0 | 4 , 5 , 6 ) = 0.00007848300723                    ( in 4 seconds with V41 )
    •  z_B( 5 , 0.8 | 4 , 5 , 6 ) = 3.052444648                             ( 1 second with V41 )
 

Notes:

-If w = n = 1 , we get the Riemann Zeta function.
-"BZF" is very slow !

-Since register P is used, don't interrupt this routine.
-Synthetic registers M N O P  may be replaced by any unused data registers.

-The following variant is shorter but also slower !
 
 
 

 
     ( 98 bytes / SIZE 2n+1 )
 
 

   Provided  R00 = n  ,   R01 = a₁ ,  ............ , Rnn = a_n

Examples:

    •  z_B( 7 , 0.8 | 4 , 5 , 6 ) = ?

    n = 3   STO 00    4  STO 01   5  STO 02   6  STO 03

     7    ENTER^
   0.8   XEQ "BZF"  >>>>  z_B( 7 , 0.8 | 4 , 5 , 6 ) = 4.768395227           ---Execution time = 94s---

-Likewise,

    •  z_B( 7 , 1.2 | 4 , 5 , 6 ) = 0.279095669                          ---Execution time = 2m57s---

Notes:

-Since register P is used, don't interrupt this routine.
-Register a is also used, so "BZF" cannot be called as more than a first level subroutine.
-But all these synthetic registers may be replaced by unused data registers...
 

     b) Integral Representation
 

-Barnes' Zeta function may also be computed by the integral:

   z_B( s,w | a₁ , ... , a_n ) = ( 1/Gamma(s) )  §₀^infinity  exp(-w.t)  t^s-1 dt / [ ( 1-exp(-a₁.t) ) ................. ( 1-exp(-a_n.t) ) ]
 

-The following program calculates the integral with Gauss-Legendre 10-point formula.
 

Data Registers:           •  R30 = n                              ( Registers R30 thru R30+nn are to be initialized before executing "BZF2" )

                                      •  R31 = a₁     •  R32 = a₂    ......................       •  R30+nn = a_n             R26 to R29: temp
 

Flags: /
Subroutines:  "GL10"  ( cf "Gaussian Integration for the HP-41" )  or another integrator... provided it doesn't disturb registers R26 , R27 , ............
                         "1/G+"  or "GAM+"  ( cf "Gamma Function for the HP-41" )

-"GL10" may be replaced by another integrator... provided it doesn't disturb registers R26 , R27 , ............
-The limits of integrations must be in registers R01 & R02
-The name of the subroutine in R00
-And the result in R04.
 
 

 
   ( 104 bytes / SIZE 31+nnn )
 
 

  Where  I_m are successive approximations of z_B( s,w | a₁ , ... , a_n )  when the interval of integration is divide in m parts.
  Provided  R30 = n  ,   R31 = a₁ ,  ............ , R30+nn = a_n

Example:    Evaluate   z_B( 5 , 0.8 | 4 , 5 , 6 )

     n = 3   STO 30    4  STO 31   5  STO 32   6  STO 33

     5    ENTER^
   0.8   XEQ "BZF2"  >>>>   I₁ = 3.225867659                          ---Execution time = 42s---

    4       R/S     >>>>    I₄ = 3.047551064
    8       R/S     >>>>    I₈ = 3.052316212
   16      R/S     >>>>   I₁₆ = 3.052446238
   32      R/S     >>>>   I₃₂ = 3.052445042
   64      R/S     >>>>   I₆₄ = 3.052445045

Notes:

-Most of these results come from V41 in turbo mode
-The last 2 estimations are certainly the best.
-Compared to the first program listed in the paragraph 1°) a)  (  z_B( 5 , 0.8 | 4 , 5 , 6 ) = 3.052444648 )
  we see that "BZF" gives an error of about -4 E-7
 

2°)  Epstein Zeta Function
 

-This function is defined by

    z_E( s,w | a₁ , ... , a_n ) = SUM_k1,....,kn  1 / ( w + a₁ k₁² + a₂ k₂² + ............... + a_n k_n² ) ^s

    with   k₁ , k₂ , .......... , k_n  =  ........ -3 , -2 , -1 , 0 , +1 , +2 , +3 , .............

    and    s > n / 2  ,  w >= 0

-If w = 0 , the term corresponding to   k₁ = k₂ = .......... = k_n = 0  is omitted
 

Data Registers:           •  R00 = n                                       ( Registers R00 thru Rnn are to be initialized before executing "EZF" )

                                      •  R01 = a₁     •  R02 = a₂    ......................       •  Rnn = a_n
                                         Rn+1 = k₁      Rn+2 = k₂  ......................          R2n = k_n

Flags: /
Subroutines: /

-Line 111 is a three-byte  GTO 01
 
 

 
        ( 173 bytes / SIZE 2n+1 )
 
 

   Provided  R00 = n  ,   R01 = a₁ ,  ............ , Rnn = a_n

Examples:

    •  z_E( 7 , 0.8 | 4 , 5 , 6 , 7 ) = ?

    n = 4  STO 00

     4  STO 01   5  STO 02   6  STO 03   7  STO 04

     7    ENTER^
    0.8  XEQ "EZF"  >>>>  4.768419976                          ---Execution time = 81s---

-Likewise, we get:

   •  z_E( 7 , 0.8 | 4 , 5 , 6 ) = 4.768418547                 in 37 seconds
   •  z_E( 7 , 1.2 | 4 , 5 , 6 ) = 0.279109441                 in 44 seconds
   •  z_E( 7 , 0 | 4 , 5 , 6 ) = 0.0001563199653           in 135 seconds

Notes:

-This program is also very slow !

-Since register P is used, don't interrupt"EZT".
-Synthetic registers M N O P  may be replaced by any unused data registers.
-Moreover, large roundoff-errors are to be expected, for example

    •  z_E( 7 , 0.8 | 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 ) = 4.803314897  ( obtained with V41 )

 whereas the correct result is about  4.803316...

-Here is also a shorter ( but slower ) version:
 
 

 
        ( 111 bytes / SIZE 2n+1 )
 
 

   Provided  R00 = n  ,   R01 = a₁ ,  ............ , Rnn = a_n

Examples:

    •  z_E( 7 , 0.8 | 4 , 5 , 6 , 7 ) = ?

    n = 4  STO 00

     4  STO 01   5  STO 02   6  STO 03   7  STO 04

     7    ENTER^
    0.8  XEQ "EZF"  >>>>  4.768419976                          ---Execution time = 122s---

Notes:

-Since register P is used, don't interrupt this routine.
-Register a is also used, so "EZF" cannot be called as more than a first level subroutine.
-But all these synthetic registers may be replaced by unused data registers...
 

3°)  Multiple Zeta Function
 

     a) Program#1
 

-According to several authors, the multiple zeta function is defined by

  z( s₁ , ... , s_n ) = SUM_{n1>n2>....,nk>0}  1 / ( n₁^s1 .......... n_k^sk )

-"MZF" uses this definition.
 

Data Registers:           •  R00 = k                                            ( Registers R00 thru Rnn are to be initialized before executing "MZF" )

                                      •  R01 = s₁     •  R02 = s₂    ......................       •  Rkk = s_k
                                         Rk+1 = n₁      Rk+2 = n₂  ......................          R2k = n_k

Flags: /
Subroutines: /
 
 

 
     ( 161 bytes / SIZE 2k+1 )
 
 

   Provided  R00 = n  ,   R01 = s₁ ,  ............ , Rkk = s_k

Examples:

    •  z( 9 , 9 , 9 ) = ?

   n = 3   STO 00          9  STO 01  STO 02  STO 03

   XEQ "MZF"  >>>>   z( 9 , 9 , 9 ) = 0.0000001081746864           ---Execution time = 109s---

-Likewise:

   •  z( 7 , 7 , 7 ) =  0.000004231442970     ( in  5mn15s )
   •  z( 7 , 8 , 9 ) =  0.000002109230767     ( in 3mn39s )
   •  z( 4 , 5 , 6 ) =  0.0006551508400         ( in a few seconds with V41 )
   •  z( 3 , 4 , 5 ) =  0.005468472467           ( in 35 seconds with V41 )

Notes:

-If k = 1 , we get the Riemann Zeta function.

-"MZF" is very slow, unless the arguments are very large. Roundoff errors are often large too !
-Since register P is used, don't interrupt the program.
-Synthetic registers M N O P  may be replaced by any unused registers
-In this case, if you replace register M by, say R11, delete line 103 and replace line 01 CLA by  CLX  STO 11
 

     b) Program#2
 

-According to other authors, the multiple zeta function is defined by

  z( s₁ , ... , s_n ) = SUM_{0<n1<n2<,...,<nk}  1 / ( n₁^s1 .......... n_k^sk )
 

-"MZF2" employs this 2nd definition.
 
 

Data Registers:           •  R00 = k                                           ( Registers R00 thru Rnn are to be initialized before executing "MZF2" )

                                      •  R01 = s₁     •  R02 = s₂    ......................       •  Rkk = s_k
                                         Rk+1 = n₁      Rk+2 = n₂  ......................          R2k = n_k

Flags: /
Subroutines: /
 
 

 
     ( 155 bytes / SIZE 2k+1 )
 
 

   Provided  R00 = n  ,   R01 = s₁ ,  ............ , Rkk = s_k

Example:   "MZF2"  returns the same results as  "MZF"  provided the s_i are stored in the reverse order, for instance:

   3   STO 00     9  STO 01   8  STO 02   7  STO 03

   XEQ "MZF2"  >>>>   z( 9 , 8 , 7 ) =  0.000002109230767     ( in 3m43s instead of 3m39s )

-Though it saves a few bytes, this second version is a little bit slower than "MZF"
-So, it seems preferable to use "MZF".
 
 

References:

[1]  Klaus Kirsten - "Basic zeta functions and some applications in physics"
[2]  I. Horozof - "Multiple Zeta Functions, Modular Forms and Adeles"
[3]  Borwein, Bradley, Broadhurst, Lisonek - "Combinatorial Aspects of Multiple Zeta Values"