hp41programs

Multidzeta

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

    zB( s,w | a1 , ... , an ) = SUMk1,....,kn=0,1,2,....    1 / ( w + a1 k1 + a2 k2 + ............... + an kn ) s

    and    s > n   ,  w >= 0
 

-If w = 0 , the term corresponding to   k1 = k2 = .......... = kn = 0  is omitted
 

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

                                      •  R01 = a1     •  R02 = a2    ......................       •  Rnn = an
                                         Rn+1 = k1      Rn+2 = k2  ......................          R2n = kn

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 = a1 ,  ............ , Rnn = an

Examples:

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

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

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

-Likewise,

    •  zB( 7 , 1.2 | 4 , 5 , 6 ) = 0.279095669                          ---Execution time = 97s---
    •  zB( 7 , 0 | 4 , 5 , 6 ) = 0.00007848300723                    ( in 4 seconds with V41 )
    •  zB( 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.
 

     b) Integral Representation
 

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

   zB( s,w | a1 , ... , an ) = ( 1/Gamma(s) )  §0infinity  exp(-w.t)  ts-1 dt / [ ( 1-exp(-a1.t) ) ................. ( 1-exp(-an.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 = a1     •  R32 = a2    ......................       •  R30+nn = an             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  Im are successive approximations of   zB( s,w | a1 , ... , an )  when the interval of integration is divide in m parts.
  Provided  R30 = n  ,   R31 = a1 ,  ............ , R30+nn = an

Example:    Evaluate   zB( 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"  >>>>   I1 = 3.225867659                          ---Execution time = 42s---

    4       R/S     >>>>    I4 = 3.047551064
    8       R/S     >>>>    I8 = 3.052316212
   16      R/S     >>>>   I16 = 3.052446238
   32      R/S     >>>>   I32 = 3.052445042
   64      R/S     >>>>   I64 = 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)  (  zB( 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

    zE( s,w | a1 , ... , an ) = SUMk1,....,kn  1 / ( w + a1 k12 + a2 k22 + ............... + an kn2 ) s

    with   k1 , k2 , .......... , kn  =  ........ -3 , -2 , -1 , 0 , +1 , +2 , +3 , .............

    and    s > n / 2  ,  w >= 0

-If w = 0 , the term corresponding to   k1 = k2 = .......... = kn = 0  is omitted
 

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

                                      •  R01 = a1     •  R02 = a2    ......................       •  Rnn = an
                                         Rn+1 = k1      Rn+2 = k2  ......................          R2n = kn

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 = a1 ,  ............ , Rnn = an

Examples:

    •  zE( 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:

   •  zE( 7 , 0.8 | 4 , 5 , 6 ) = 4.768418547                 in 37 seconds
   •  zE( 7 , 1.2 | 4 , 5 , 6 ) = 0.279109441                 in 44 seconds
   •  zE( 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

    •  zE( 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...
 

3°)  Multiple Zeta Function
 

     a) Program#1
 

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

  z( s1 , ... , sn ) = SUMn1>n2>....,nk>0  1 / ( n1s1 .......... nksk )

-"MZF" uses this definition.
 

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

                                      •  R01 = s1     •  R02 = s2    ......................       •  Rkk = sk
                                         Rk+1 = n1      Rk+2 = n2  ......................          R2k = nk

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 = s1 ,  ............ , Rkk = sk

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( s1 , ... , sn ) = SUM0<n1<n2<,...,<nk  1 / ( n1s1 .......... nksk )
 

-"MZF2" employs this 2nd definition.
 
 

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

                                      •  R01 = s1     •  R02 = s2    ......................       •  Rkk = sk
                                         Rk+1 = n1      Rk+2 = n2  ......................          R2k = nk

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 = s1 ,  ............ , Rkk = sk

Example:   "MZF2"  returns the same results as  "MZF"  provided the si 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"