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: /
 
 

01  LBL "BZF"    02  STO N  03  X<>Y  04  CHS  05  STO O  06  CHS  07  Y^X  08  X#0?  09  1/X  10  STO M  11  RCL 00   12  ST+ X  13   E3  14  /  15  RCL 00  16  +  17  ISG X  18  CLRGX

19  RCL 00  20  .1  21  %  22  +  23  RCL 00  24  +  25  STO P  26  RCL N  27  RDN  28  GTO 03  29  LBL 01  30  RCL P  31  RCL 00   32  -  33  RCL IND X  34  ST+ IND P  35  RCL N  36  RCL IND P

37  +  38  STO Z  39  RCL O  40  Y^X  41  RCL M  42  +  43  STO M  44  LASTX  45  X=Y?  46  GTO 04  47  R^  48  RCL P  49  FRC  50  RCL X  51  ST+ Y  52  RCL 00        53  ST+ X  54  +

55  X<> P  56  INT  57  +  58  RCL IND X  59  LBL 02  60  STO IND Y  61  ISG Y  62  GTO 02   63  R^  64  LBL 03  65  R^  66  RCL IND 00  67  +  68  STO Z  69  RCL O   70  Y^X  71  RCL M  72  +

73  STO M  74  LASTX  75  X#Y?  76  GTO 03   77  LBL 04  78  DSE P  79  GTO 01  80  X<> N  81  SIGN  82  X<> O           83  CHS  84  X<>Y  85  CLA  86  END

 
     ( 134 bytes / SIZE 2n+1 )
 
 

STACK	INPUTS	OUTPUTS
Y	s	s
X	w	zB( s,w | a1 , ... , an )
L	/	w

   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 !
 
 
 

01 LBL "BZF"     02 STO N  03 X<>Y  04 CHS  05 STO O  06 CHS  07 Y^X  08 X#0?  09 1/X  10 STO M  11 RCL 00  12  E3

13 /  14 STO a  15 ST+ X  16 RCL 00  17 ST+ a  18 ST+ a  19 +  20 ISG X  21 CLRGX  22 LBL 01          23 RCL a  24 STO P

25 LBL 02  26 CLX  27 SIGN  28 ST+ IND P  29 RCL 00  30 ST+ X  31 RCL N  32 LBL 03  33 RCL IND Y  34 RCL IND 00  35 *  36 +

37 DSE Y  38 DSE 00  39 GTO 03  40 X<>Y  41 STO 00  42 CLX  43 RCL O  44 Y^X  45 RCL M          46 +  47 STO M  48 LASTX

49 X#Y?  50 GTO 01  51 CLX  52 STO IND P   53 DSE P  54 GTO 02  55 X<> N  56 SIGN  57 X<> O  58 CHS  59 X<>Y  60 CLA  61 END

 
     ( 98 bytes / SIZE 2n+1 )
 
 

STACK	INPUTS	OUTPUTS
Y	s	s
X	w	zB( s,w | a1 , ... , an )
L	/	w

   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.
 
 

01  LBL "BZF2"   02  DEG  03  STO 28  04  X<>Y  05  STO 29  06  "T"  07  ASTO 00  08  CLX  09  STO 01   10  90  11  STO 02  12  30.03

13  RCL 30  14  +  15  STO 26  16  SIGN  17  LBL 01  18  STO 03  19  XEQ "GL10"  20  RCL 29  21  XEQ "1/G+"  22  RCL 04   23  ABS  24  *

25  D-R  26  RTN  27  GTO 01  28  LBL "T"  29  TAN  30  STO Y  31  RCL 29         32  1  33  -  34  Y^X  35  X<>Y  36  X^2

37  ENTER^  38  SIGN  39  +  40  *  41  X<>Y  42  RCL 28   43  *  44  CHS  45  E^X  46  *  47  RCL 26          48  STO 27

49  RDN  50  LBL 02  51  X<>Y  52  RCL IND 27  53  *  54  CHS  55  E^X-1  56  /  57  DSE 27  58  GTO 02   59  RTN  60  END

 
   ( 104 bytes / SIZE 31+nnn )
 
 

STACK	INPUTS	OUTPUT1	INPUT2	OUTPUT2
Y	s	/	/	/
X	w	I1	m	Im

  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
 
 

01  LBL "EZF"     02  STO N   03  X<>Y   04  STO O   05  Y^X   06  X#0?   07  1/X   08  STO M   09  RCL 00    10  ST+ X   11   E3   12  /   13  RCL 00   14  +   15  ISG X   16  CLRGX   17  2   18  RCL 00   19  .1   20  %   21  +   22  RCL 00   23  +   24  STO P

25  RCL N   26  X<>Y   27  ENTER^   28  GTO 03    29  LBL 01   30  R^   31  RCL IND P    32  X#0?   33  SIGN   34  ENTER^   35  SIGN   36  ST+ IND P   37  +   38  /   39  STO Q   40  RCL P   41  INT   42  STO Y   43  RCL 00   44  -   45  STO 00   46  CLX   47  RCL N   48  LBL 02

49  RCL IND Y   50  X^2   51  RCL IND 00   52  *   53  +   54  DSE Y   55  DSE 00   56  GTO 02   57  X<>Y   58  STO 00    59  X<> Q   60  X<>Y   61  RCL X   62  RCL O   63  Y^X   64  R^   65  X<>Y   66  /   67  RCL M   68  +   69  STO M   70  LASTX   71  X=Y?   72  GTO 04

73  RDN   74  X<> Z   75  ST+ X   76  STO Z   77  CLX   78  RCL P   79  FRC   80  RCL 00    81  ST+ X   82  +   83  STO P   84  ENTER^   85  LBL 03   86  SIGN   87  ST+ IND P   88  RDN   89  CLX   90  RCL IND P   91  X^2   92  RCL IND 00   93  *   94  RCL Y   95  +   96  RCL O

97  Y^X   98  R^   99  X<>Y 100  / 101  RCL M 102  + 103  STO M 104  LASTX 105  X#Y? 106  GTO 03  107  LBL 04 108  CLX 109  STO IND P  110  DSE P 111  GTO 01 112  X<> N 113  SIGN 114  X<> O 115  X<>Y 116  CLA 117  END

 
        ( 173 bytes / SIZE 2n+1 )
 
 

STACK	INPUTS	OUTPUTS
Y	s	s
X	w	zE( s,w | a1 , ... , an )
L	/	w

   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:
 
 

01 LBL "EZF"     02 STO N  03 X<>Y  04 STO O  05 Y^X  06 X#0?  07 1/X  08 STO M  09 RCL 00  10  E3  11 /  12 STO a  13 ST+ X  14 RCL 00

15 ST+ a  16 ST+ a  17 +  18 ISG X  19 CLRGX  20 LBL 01  21 RCL a  22 STO P  23 LBL 02  24 CLX  25 SIGN  26 ST+ IND P   27 STO Q  28 RCL 00

29 ST+ X  30 RCL N  31 LBL 03  32 RCL IND Y  33 X=0?  34 GTO 03  35 X<> Q  36 ST+ X  37 X<> Q  38 X^2  39 LBL 03  40 RCL IND 00  41 *  42 +

43 DSE Y  44 DSE 00  45 GTO 03  46 X<>Y  47 STO 00        48 X<> Q  49 X<>Y  50 RCL O  51 Y^X  52 /  53 RCL M  54 +  55 STO M

56 LASTX  57 X#Y?  58 GTO 01  59 CLX  60 STO IND P   61 DSE P  62 GTO 02  63 X<> N  64 SIGN  65 X<> O  66 X<>Y  67 CLA  68 END

 
        ( 111 bytes / SIZE 2n+1 )
 
 

STACK	INPUTS	OUTPUTS
Y	s	s
X	w	zE( s,w | a1 , ... , an )
L	/	w

   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: /
 
 

01  LBL "MZF"   02  CLA   03  RCL 00   04  ENTER^   05  ST+ X   06   E3   07  /   08  +   09  ISG X   10  STO N   11  STO O   12  INT   13  RCL 00    14  LBL 00   15  STO IND Y   16  ISG Y   17  CLX   18  DSE X   19  GTO 00   20  SIGN   21  ST- IND N

22  2   23  STO P   24  RCL N   25  ENTER^   26  SIGN   27  LBL 01   28  ISG Y   29  X<0?   30  GTO 02   31  RCL IND Y   32  RCL IND P   33  Y^X   34  /   35  ISG P   36  CLX   37  GTO 01    38  LBL 02   39  RDN   40  GTO 06   41  LBL 03   42  RCL N

43  INT   44  RCL 00   45   E3   46  /   47  +   48  RCL IND N   49  LBL 04   50  ENTER^   51  SIGN   52  +   53  STO IND Y   54  DSE Y   55  GTO 04    56  SIGN   57  ST+ X   58  STO P   59  RCL O   60  ENTER^   61  SIGN   62  ST+ Y   63  LBL 05

64  RCL IND Y   65  RCL IND P   66  Y^X   67  /   68  ISG P   69  CLX   70  ISG Y   71  GTO 05    72  STO Y   73  RCL IND O   74  RCL 01   75  Y^X   76  /   77  RCL M   78  +   79  STO M   80  LASTX   81  X=Y?   82  GTO 07   83  RCL O   84  STO N

85  LBL 06   86  SIGN   87  ST+ IND N   88  R^   89  STO Y   90  RCL IND N   91  RCL 01    92  Y^X   93  /   94  RCL M   95  +   96  STO M   97  LASTX   98  X#Y?   99  GTO 06 100  LBL 07 101  ISG N 102  GTO 03 103  CLA 104  END

 
     ( 161 bytes / SIZE 2k+1 )
 
 

STACK	INPUTS	OUTPUTS
X	/	z( s1 , ... , sk )

   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: /
 
 

01  LBL "MZF2"  02  CLA  03  RCL 00  04  ENTER^  05  ST+ Y  06   E3  07  /  08  +  09  STO N  10  STO O  11  INT  12  RCL 00   13  LBL 00   14  STO IND Y  15  DSE Y  16  DSE X  17  GTO 00  18  SIGN  19  ST- IND N  20  RCL 00

21  1  22  LBL 01  23  DSE Y  24  X<0?  25  GTO 02  26  RCL Y  27  RCL IND Z  28  Y^X  29  /  30  GTO 01   31  LBL 02  32  RDN  33  GTO 06  34  LBL 03  35  RCL O  36  INT  37   E3  38  /  39  RCL N  40  INT

41  +  42  RCL IND X  43  LBL 04  44  ENTER^  45  SIGN  46  +  47  STO IND Y  48  ISG Y  49  GTO 04  50  RCL 00   51  STO P  52  RCL O  53  INT  54  ENTER^  55  SIGN  56  ST- Y  57  ST- P  58  LBL 05  59  RCL IND Y  60  RCL IND P

61  Y^X  62  /  63  DSE Y  64  DSE P  65  GTO 05  66  STO Y  67  RCL IND O  68  RCL IND 00  69  Y^X  70  /  71  RCL M  72  +  73  STO M  74  LASTX  75  X=Y?  76  GTO 07   77  RCL O  78  STO N  79  LBL 06  80  SIGN

81  ST+ IND N  82  R^  83  STO Y  84  RCL IND N  85  RCL IND 00  86  Y^X  87  /  88  RCL M  89  +  90  STO M  91  LASTX  92  X#Y?  93  GTO 06   94  LBL 07  95  DSE N  96  GTO 03  97  CLA  98  END

 
     ( 155 bytes / SIZE 2k+1 )
 
 

STACK	INPUTS	OUTPUTS
X	/	z( s1 , ... , sk )

   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"