Riemann Zeta Function for the HP-41

Overview
 

 1°) Real Arguments

  a)  Elementary Method
  b)  Borwein Algorithm

 2°) Complex Arguments
 3°) Quaternion Arguments

  a)  Elementary Method
  b)  Borwein Algorithm
 

1°)  Real Arguments
 

     a) Elementary Method
 

-The Riemann Zeta function is defined by the series:       Zeta(x) = 1 + 1/2^x + 1/3^x + ........ + 1/n^x + ......     ( for x > 1)
  but the following program uses the formula:          Zeta(x) = ( 1 - 1/2^x + 1/3^x - 1/4^x + ...... ) / ( 1 - 2^1-x )   where roundoff errors are smaller.

-If  x < 0.5 , we have:   Zeta(x) = GAM(1-x)  2^x  PI^x-1  Sin(PI.x/2)  Zeta(1-x)
 
 

Data Registers:  R00 thru R03 ( temp )
Flags: /
Subroutine:  "GAM"  ( cf "Gamma Function for the HP-41" )  if x < 1 only.
 
 
 

01  LBL "ZETA" 02  1 03  X<>Y 04  X>Y? 05  GTO 00       06  STO 02 07  - 08  STO 03 09  XEQ 00 10  1 11  ASIN 12  RCL 02

13  * 14  SIN 15  * 16  PI 17  RCL 03       18  Y^X 19  / 20  2 21  RCL 02 22  Y^X 23  * 24  X<> 03

25  XEQ "GAM" 26  RCL 03        27  * 28  RTN 29  LBL 00 30  CHS 31  STO 01  32  1 33  STO 00 34  ENTER^ 35  LBL 01 36  CLX

37  SIGN 38  RCL 00 39  + 40  STO 00        41  RCL 01  42  Y^X 43  - 44  CHS 45  LASTX 46  RND  47  X#0? 48  GTO 01

49  SIGN 50  2 51  RCL 01        52  Y^X 53  ST+ X 54  - 55  / 56  ABS 57  END

 
   ( 76 bytes / SIZE 004 )
 
 

STACK	INPUT	OUTPUT
X	x	Zeta(x)

 
Examples:

  FIX 7       3    XEQ "ZETA"    returns    Zeta(3) = 1.2020569          ( in 3m14s )
  FIX 9       3           R/S            -------    Zeta(3) = 1.202056903      ( in 15m12s )

  FIX 9   -7.49        R/S            -------  Zeta(-7.49) = 0.003312040169  ( in 15s )

Notes:

-The accuracy is determined by the display format ( line 46 ).
-Execution time and Zeta(x) tend to infinity as x tends to 1.
 

     b) Borwein Algorithm
 

-This second version uses the method described in reference [3]
 

Data Registers:  R00 to R07: temp
Flags: /
Subroutine:    "1/G+"  ( cf "Gamma function for the HP-41" )

-Line 69 is a synthetic TEXT0 instruction ( NOP )
 
 

01  LBL "ZETA"  02  X#0?  03  GTO 00        04  .5  05  CHS  06  RTN  07  LBL 00  08  STO 01  09  STO 07  10  .5  11  X<=Y?  12  GTO 00  13  SIGN  14  X<>Y  15  -  16  STO 01  17  XEQ 00  18  1

19  ASIN  20  RCL 07        21  *  22  SIN  23  *  24  PI  25  ST/ Y  26  ST+ X  27  RCL 07   28  Y^X  29  *  30  X<> 01  31  XEQ "1/G+"  32  ST/ 01  33  X<> 01  34  RTN  35  LBL 00  36  14

37  STO 00  38  STO 03  39  SIGN  40  STO 02        41  STO 04  42  STO 05  43  CHS  44  STO 06   45  CLX  46  LBL 01  47  RCL 05  48  RCL 00  49  RCL 01  50  CHS  51  Y^X  52  *  53  RCL 06  54  CHS

55  STO 06  56  *  57  +  58  RCL 00        59  ENTER^  60  ST+ Y  61  ST* Y  62  -  63  RCL 04   64  *  65  RCL 03  66  X^2  67  RCL 00  68  DSE X  69  ""  70  X^2  71  -  72  ST+ X

73  /  74  STO 04        75  ST+ 05  76  X<>Y  77  DSE 00  78  GTO 01  79  RCL 05   80  /  81  2  82  LN  83  1  84  RCL 01  85  -  86  *  87  E^X-1  88  /  89  END

 
    ( 122 bytes / SIZE 008 )
 
 

STACK	INPUT	OUTPUT
X	x	Zeta(x)

 
Examples:

      3    XEQ "ZETA"   >>>>      Zeta(3)    = 1.202056903                   ---Execution time = 21s---
  -7.49 XEQ "ZETA"   >>>>    Zeta(-7.49) = 0.003312040169            ---Execution time = 24s---
    1.1   XEQ "ZETA"  >>>>      Zeta(1.1)  = 10.58444847                   ---Execution time = 21s---

Notes:

  Zeta(0) = -1/2
-Borwein Algorithm is really fantastic !
-You get  Zeta(1.001) in 21 seconds, which seemed impossible before...
 

2°)  Complex Arguments
 

-"ZETAZ" also uses Borwein algorithm if Re(z) > 0.5

-If  Re(z) < 0.5 , we've used    Zeta(z) = Zeta(1-z)  Pi ^z-1/2  Gamma((1-z)/2) / Gamma(z/2)

-The other formula may also be used but it seems to produce more roundoff-errors.
 

Data Registers:  R00 to R11: temp      When the program stops, the result is in R04-R05
Flags: /
Subroutines:   "GAMZ+"  ( cf "Gamma Function for the HP-41" )
                           "Z*Z"  "Z/Z"  "Z^Z"  ( cf "Complex Functions for the HP-41" )
 

-Line 128 is a synthetic TEXT0 instruction ( NOP )
 
 
 

01  LBL "ZETAZ"   02  X=0?   03  X#Y?   04  GTO 00               05  .5   06  CHS   07  RTN   08  LBL 00   09  X<>Y   10  STO 07    11  STO 11   12  X<>Y   13  STO 06   14  STO 10   15  .5   16  X<=Y?   17  GTO 00   18  SIGN   19  X<>Y   20  -   21  STO 06   22  RCL 07   23  CHS   24  STO 07   25  XEQ 00   26  RCL 11   27  CHS   28  1   29  RCL 10   30  -

31  2   32  ST/ Z   33  /   34  XEQ "GAMZ+"   35  RCL 05               36  RCL 04   37  XEQ "Z*Z"   38  STO 04   39  X<>Y   40  STO 05   41  PI   42  RCL 11   43  RCL 10   44  .5   45  ST* 10   46  ST* 11   47  -   48  0   49  RDN   50  XEQ "Z^Z"   51  RCL 05   52  RCL 04   53  XEQ "Z*Z"   54  STO 04   55  X<>Y   56  STO 05   57  RCL 11   58  RCL 10   59  XEQ "GAMZ+"   60  RCL 05

61  RCL 04   62  GTO 02   63  LBL 00   64  PI   65  2   66  /   67  RCL 07               68  ABS   69  ST* Y   70  ST+ X   71  LN1+X   72  +   73  3 E10   74  LN   75  +   76  8   77  SQRT   78  3   79  +   80  LN   81  /   82  INT   83  1   84  +   85  STO 00   86  STO 02   87  LASTX   88  STO 01   89  STO 03   90  STO 04

91  CHS   92  X<>Y   93  Y^X   94  CHS   95  STO 05   96  CLX   97  STO 08               98  STO 09   99  LBL 01 100  CLST 101  RCL 00 102  RCL 07 103  CHS 104  RCL 06 105  CHS 106  XEQ "Z^Z" 107  RCL 04 108  RCL 05 109  CHS 110  STO 05 111  * 112  ST* Z 113  * 114  ST+ 08 115  X<>Y 116  ST+ 09 117  RCL 00 118  ENTER^ 119  ST+ Y 120  ST* Y

121  - 122  RCL 03 123  * 124  RCL 02 125  X^2 126  RCL 00             127  DSE X 128  "" 129  X^2 130  - 131  ST+ X 132  / 133  STO 03 134  ST+ 04 135  DSE 00 136  GTO 01 137  RCL 04 138  ST/ 08 139  ST/ 09 140  2 141  LN 142  1 143  RCL 06 144  - 145  * 146  E^X-1 147  STO 01 148  RCL 07 149  CHS 150  2

151  LN 152  * 153  1 154  RAD 155  P-R 156  ENTER^ 157  DEG 158  1 159  - 160  RCL 01             161  ST* Z 162  X<> T 163  ST* T 164  ST+ T 165  RDN 166  + 167  RCL 09 168  RCL 08 169  LBL 02 170  R^ 171  R^ 172  XEQ "Z/Z" 173  STO 04 174  X<>Y 175  STO 05 176  X<>Y 177  END

 
    ( 251 bytes / SIZE 012 )
 
 

STACK	INPUTS	OUTPUTS
Y	y	y'
X	x	x'

       Zeta(x+i.y) = x' + i.y'

Example:

   7   ENTER^
  -6  XEQ "ZETAZ"  >>>>  -4.135010534
                                RDN    -1.198478879

-So,  Zeta(-6+7.i) = -4.135010534 - 1.198478879 i
 

3°)  Quaternion Arguments
 

     a) Elementary Method
 

-"ZETAQ" uses  Zeta(q) = 1 + 1/2^q + 1/3^q + ........ + 1/n^q + ......     ( for Re(q) > 1)

-And if Re(q) < 0 ,    Zeta(q) = Zeta(1-q)  Pi ^q-1/2  Gamma((1-q)/2) / Gamma(q/2)
 

Data Registers:  R00 to R24: temp
Flags: /
Subroutines:     "GAMQ"  ( cf "Quaternionic Special Functions for the HP-41" )
                            "Q*Q"  "1/Q"  "E^Q"  ( cf "Quaternions for the HP-41" )

-Line 16 is a three-byte GTO 00
 
 
 

01  LBL "ZETAQ"   02  STO 05              03  STO 17   04  RDN   05  STO 06   06  STO 18    07  RDN   08  STO 07   09  STO 19   10  X<>Y   11  STO 08   12  STO 20   13  RCL 05   14  1   15  X<Y?   16  GTO 00   17  CHS   18  ST* 05   19  ST* 06   20  ST* 07   21  ST* 08   22  ST- 05   23  XEQ 00   24  STO 21   25  RDN   26  STO 22   27  RDN   28  STO 23   29  X<>Y   30  STO 24

31  2   32  ST/ 05   33  ST/ 06   34  ST/ 07   35  ST/ 08   36  RCL 08   37  RCL 07   38  RCL 06   39  RCL 05   40  XEQ "GAMQ"   41  STO 05              42  RDN   43  STO 06   44  RDN   45  STO 07   46  X<>Y   47  STO 08   48  21.001004   49  REGMOVE   50  XEQ "Q*Q"   51  STO 21   52  RDN   53  STO 22   54  RDN   55  STO 23   56  X<>Y   57  STO 24   58  RCL 17   59  SIGN   60  RCL 20

61  RCL 19   62  RCL 18   63  2   64  ST/ L   65  ST/ Y   66  ST/ Z   67  ST/ T   68  X<> L   69  XEQ "GAMQ"   70  XEQ "1/Q"   71  STO 05              72  RDN   73  STO 06    74  RDN   75  STO 07   76  X<>Y   77  STO 08   78  21.001004   79  REGMOVE   80  XEQ "Q*Q"   81  STO 05   82  RDN   83  STO 06   84  RDN   85  STO 07   86  X<>Y   87  STO 08   88  .5   89  ST- 17   90  RCL 20

91  RCL 19   92  RCL 18   93  PI   94  LN   95  ST* 17   96  ST* T   97  ST* Z   98  *   99  RCL 17 100  XEQ "E^Q" 101  STO 01            102  RDN 103  STO 02  104  RDN 105  STO 03 106  X<>Y 107  STO 04 108  XEQ "Q*Q" 109  RTN 110  LBL 00 111  CLX 112  STO 02 113  STO 03 114  STO 04 115  SIGN 116  STO 00 117  STO 01 118  LBL 01 119  ISG 00 120  CLX

121  RCL 08 122  RCL 07 123  RCL 06 124  RCL 00 125  LN 126  CHS 127  SIGN 128  CLX 129  RCL 05            130  X<> L 131  ST* L 132  ST* Y 133  ST* Z 134  ST* T 135  X<> L 136  XEQ "E^Q" 137  STO 09 138  CLX 139  RCL 02 140  + 141  STO 02 142  LASTX 143  - 144  ABS 145  X<>Y 146  RCL 03 147  + 148  STO 03 149  LASTX 150  -

151  ABS 152  + 153  X<>Y 154  RCL 04 155  + 156  STO 04            157  LASTX 158  - 159  ABS 160  + 161  RCL 09 162  RCL 01 163  + 164  STO 01 165  LASTX 166  - 167  ABS 168  + 169  X#0? 170  GTO 01 171  RCL 04 172  RCL 03 173  RCL 02 174  RCL 01 175  END

 
   ( 288 bytes / SIZE 025 )
 
 

STACK	INPUTS	OUTPUTS
T	t	t'
Z	z	z'
Y	y	y'
X	x	x'

       Zeta(x+i.y+j.z+k.t) = x' + i.y' + j.z' + k.t'

Example:

   7.9   ENTER^
   7.8   ENTER^
   7.7   ENTER^
 -7.6   XEQ "ZETAQ"    >>>>     762.9852911                                    ---Execution time = 2mn---
                                      RDN      -14.26823437
                                      RDN      -14.45353627
                                      RDN      -14.63883820

-Thus,   Zeta ( -7.6 + 7.7 i + 7.8 j + 7.9 k ) = 762.9852911 - 14.26823437 i - 14.45353627 j - 14.63883820 k

Notes:

-Lines 137 to 168 may be replaced by  DSQ   1    ( cf "M-Code routines for hypercomplex numbers" )
-The execution time corresponds to this modification...

-Add  RND  after line 168 and the accuracy will be controlled by the display settings.
 

     b) Borwein Algorithm
 
 

Data Registers:  R00 to R28: temp
Flags: /
Subroutines:     "GAMQ"  ( cf "Quaternionic Special Functions for the HP-41" )
                            "Q*Q"  "1/Q"  "E^Q"  ( cf "Quaternions for the HP-41" )

-Lines 30 and 137 are three-byte GTOs
 
 

01  LBL "ZETAQ"   02  STO 05              03  STO 17   04  X^2   05  X<>Y   06  STO 06   07  STO 18   08  X^2   09  +   10  X<>Y   11  STO 07   12  STO 19   13  X^2   14  +   15  X<>Y   16  STO 08   17  STO 20   18  X^2   19  +   20  X#0?   21  GTO 00   22  CLST   23  .5   24  CHS   25  RTN   26  LBL 00   27  RCL 05   28  .5   29  X<=Y?   30  GTO 00   31  SIGN   32  CHS   33  ST* 05   34  ST* 06   35  ST* 07   36  ST* 08   37  ST- 05   38  RCL 05   39  STO 25   40  RCL 06   41  STO 26   42  RCL 07   43  STO 27   44  RCL 08   45  STO 28   46  XEQ 00   47  STO 21   48  RDN

49  STO 22   50  RDN   51  STO 23   52  X<>Y   53  STO 24   54  RCL 28              55  RCL 27   56  RCL 26   57  2   58  ST/ 25   59  ST/ T   60  ST/ Z   61  /   62  RCL 25   63  XEQ "GAMQ"   64  STO 05   65  RDN   66  STO 06   67  RDN   68  STO 07   69  X<>Y   70  STO 08   71  RCL 21   72  STO 01   73  RCL 22   74  STO 02   75  RCL 23   76  STO 03   77  RCL 24   78  STO 04   79  XEQ "Q*Q"   80  STO 21   81  RDN   82  STO 22   83  RDN   84  STO 23   85  X<>Y   86  STO 24   87  RCL 20   88  RCL 19   89  RCL 18   90  RCL 17   91  SIGN   92  CLX   93  2   94  ST/ L   95  ST/ Y   96  ST/ Z

97  ST/ T   98  X<> L   99  XEQ "GAMQ" 100  XEQ "1/Q" 101  STO 05            102  RDN 103  STO 06 104  RDN 105  STO 07 106  X<>Y 107  STO 08 108  RCL 21 109  STO 01 110  RCL 22 111  STO 02 112  RCL 23 113  STO 03 114  RCL 24 115  STO 04 116  XEQ "Q*Q" 117  STO 01 118  RDN 119  STO 02 120  RDN 121  STO 03 122  X<>Y 123  STO 04 124  .5 125  ST- 17 126  PI 127  LN 128  ST* 17 129  ST* 18 130  ST* 19 131  ST* 20 132  RCL 20 133  RCL 19 134  RCL 18 135  RCL 17 136  XEQ "E^Q" 137  GTO 02 138  LBL 00 139  PI 140  2 141  / 142  RCL 06 143  X^2 144  RCL 07

145  X^2 146  RCL 08 147  X^2 148  + 149  + 150  SQRT 151  ST* Y 152  ST+ X 153  LN1+X 154  + 155  3 E10 156  LN 157  + 158  8 159  SQRT 160  3 161  + 162  LN 163  / 164  INT 165  1 166  + 167  STO 00            168  STO 02 169  LASTX 170  STO 01 171  STO 03 172  STO 04 173  CHS 174  X<>Y 175  Y^X 176  CHS 177  STO 13 178  CLX 179  STO 09 180  STO 10 181  STO 11 182  STO 12 183  LBL 01 184  RCL 00 185  LN 186  SIGN 187  RCL 08 188  CHS 189  RCL 07 190  CHS 191  RCL 06 192  CHS

193  RCL 05 194  CHS 195  X<> L 196  ST* L 197  ST* Y 198  ST* Z 199  ST* T 200  X<> L 201  XEQ "E^Q" 202  STO 14 203  RDN 204  STO 15            205  CLX 206  RCL 04 207  RCL 13 208  CHS 209  STO 13 210  * 211  ST* 14 212  ST* 15 213  ST* Z 214  * 215  ST+ 11 216  X<>Y 217  ST+ 12 218  RCL 14 219  ST+ 09 220  RCL 15 221  ST+ 10 222  RCL 00 223  ENTER^ 224  ST+ Y 225  ST* Y 226  - 227  RCL 03 228  * 229  RCL 02 230  X^2 231  RCL 00 232  1 233  - 234  X^2 235  - 236  ST+ X 237  / 238  STO 03 239  ST+ 04 240  DSE 00

241  GTO 01 242  RCL 09 243  RCL 04 244  / 245  STO 01 246  RCL 10 247  LASTX 248  / 249  STO 02 250  RCL 11            251  LASTX 252  / 253  STO 03 254  RCL 12 255  LASTX 256  / 257  STO 04 258  RCL 08 259  RCL 07 260  RCL 06 261  2 262  LN 263  CHS 264  ST* 05 265  ST- 05 266  ST* T 267  ST* Z 268  * 269  RCL 05 270  XEQ "E^Q" 271  SIGN 272  ST* X 273  ST- L 274  X<> L 275  XEQ "1/Q" 276  LBL 02 277  STO 05 278  RDN 279  STO 06 280  RDN 281  STO 07 282  X<>Y 283  STO 08 284  XEQ "Q*Q" 285  END

 
    ( 432 bytes / SIZE 029 )
 
 

STACK	INPUTS	OUTPUTS
T	t	t'
Z	z	z'
Y	y	y'
X	x	x'

       Zeta(x+i.y+j.z+k.t) = x' + i.y' + j.z' + k.t'

Examples:

   7.9   ENTER^
   7.8   ENTER^
   7.7   ENTER^
 -7.6   XEQ "ZETAQ"    >>>>     762.9852917                          ---Execution time = 2m45s---
                                      RDN      -14.26823447
                                      RDN      -14.45353622
                                      RDN      -14.63883832

           Zeta ( -7.6 + 7.7 i + 7.8 j + 7.9 k ) = 762.9852917 - 14.26823447 i - 14.45353622 j - 14.63883832 k
 

-Likewise,  Zeta( 1.1 + 7 i + 8 j + 9 k ) = 0.389296030 - 0.027737071 i - 0.031699509 j - 0.035661948 k          ( in 2m06s )
 

Notes:

-If q is very close to 1, there could be a loss of accuracy when 2^(1-q)-1 is executed:
-Replace lines 258 to 275 by

   1                   STO 09          RCL 16        ST* Y             RCL 10          CHS               ST* 06           RCL 06
   RCL 05         *                    *                   1                     E^X               RCL 16           ST* 07          RCL 05
   -                    STO 10          STO 09        -                     RCL 09          X=0?               ST* 08          DEG
   2                    E^X-1            RAD            +                    SIN                SIGN              RCL 08          XEQ "1/Q"
   LN                RCL 09          COS            STO 05           *                     /                     RCL 07

  and add   STO 16 after line 150.

-The last decimals may depend on the version of "E^Q" , "Q*Q" , ....  that you are using.
-If Re(q) is large, it could be preferable to employ the 1st version of "ZETAQ".
-In all other cases, Borwein algorithm is unbeatable !
 

References:

[1]  Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]   http://functions.wolfram.com/
[3]  P. Borwein - "An Efficient Algorithm for the Riemann Zeta Function"