Riemann's Zeta Function - another approach (RPL)

[Dieter]

(06-30-2017 11:06 AM)Gerson W. Barbosa Wrote:  I tested one by one with ever higher n until I got all 10 significant digits right.

LOL – I did the same with Excel (VBA) and 15-digit precision, with a threshold of 0,5 ULP compared to the (more or less) exact result with 200 terms. I now have slightly modified the estimate for the required number of terms, see below.

(06-30-2017 11:06 AM)Gerson W. Barbosa Wrote:  Do you intend to write a full range HP-41 version? I'll be traveling the next few days and I won't take any 41 along.

It seems that this would require some kind of Gamma routine, yes?

(06-30-2017 11:06 AM)Gerson W. Barbosa Wrote:  Thanks again for your great contributions, especially for your striving to always get as many significant digits as possible.

While we're at it: I just looked at JMB's routines and found two interesting things: first, in the Borwein program he uses the same e^x–1 method as my own program, and the other version has is a nice trick for evaluating series' with alternating signs.

Finally I now got a good approximation for x between 1 and 1,05 (!) which (sufficient precision provided) is within ±0,2 ULP of 10 digits:

Zeta(x) ~ 1/u + u/(0,9246 · u + 13,733) + 0,577215654
where u = x–1

All this leads to the following new and improved version:

Code:

01  LBL"ZETA"
02  STO 00
03  1
04  -
05  LN
06  ,05
07  LastX
08  x>y?
09  GTO 00
10  1/x
11  LastX
12  LastX
13  ,9246
14  *
15  13,733
16  +
17  /
18  ,577215654
19  +
20  +
21  GTO 02
22  LBL 00
23  25
24  RCL 00
25  /
26  INT
27  2
28  +
29  ST+ X
30  STO 01
31  RCL 00
32  CHS
33  STO 00
34  CLX
35  LBL 01
36  RCL Y
37  RCL 00
38  y^x
39  -
40  CHS
41  DSE Y
42  GTO 01
43  RCL 00
44  ST+ X
45  1
46  -
47  24
48  /
49  RCL 01
50  x^2
51  /
52  1
53  RCL 00
54  -
55  8
56  /
57  RCL 01
58  /
59  +
60  ,5
61  +
62  RCL 01
63  +
64  RCL 00
65  y^x
66  2
67  /
68  +
69  RCL 00
70  1
71  +
72  2
73  LN
74  *
75  e^x-1
76  CHS
77  /
78  LBL 02
79  END

Again, this works for x > 1.

Addendum for the 15C and others: The benefit of the e^x–1 method for the resulting accuracy can also be achieved in another way.
2^x/(2^x–2) – 1 = 1/(2^x–1–1)

Edit: not sure if this is a good idea for x close to 1 – see further below for an e^x–1 implementation on the 15C.

So the final multiplication with this power-of-2-fraction can also be coded like this (assuming –x is stored in R1, as in your 15C program):

Code:

...
2
RCL 1
CHS
1
-
y^x
1
-
/
+

Edit: Finally, here is a 15C version. Please also note the alternative solution with an e^x–1 implementation a few posts below.

Code:

# --------------------------------------------
# HEWLETT·PACKARD 15C Simulator program
# Created with version 3.2.00
# --------------------------------------------
# --------------------------------------------

   000 {             } 
   001 {    42 21 11 } f LBL A
   002 {       44  0 } STO 0
   003 {           1 } 1
   004 {          30 } −
   005 {       43 12 } g LN
   006 {          48 } .
   007 {           0 } 0
   008 {           5 } 5
   009 {       43 36 } g LSTx
   010 {    43 30  7 } g TEST x>y
   011 {       22  0 } GTO 0
   012 {          15 } 1/x
   013 {       43 36 } g LSTx
   014 {       43 36 } g LSTx
   015 {          48 } .
   016 {           9 } 9
   017 {           2 } 2
   018 {           4 } 4
   019 {           6 } 6
   020 {          20 } ×
   021 {           1 } 1
   022 {           3 } 3
   023 {          48 } .
   024 {           7 } 7
   025 {           3 } 3
   026 {           3 } 3
   027 {          40 } +
   028 {          10 } ÷
   029 {          48 } .
   030 {           5 } 5
   031 {           7 } 7
   032 {           7 } 7
   033 {           2 } 2
   034 {           1 } 1
   035 {           5 } 5
   036 {           6 } 6
   037 {           5 } 5
   038 {           4 } 4
   039 {          40 } +
   040 {          40 } +
   041 {       43 32 } g RTN
   042 {    42 21  0 } f LBL 0
   043 {           2 } 2
   044 {           5 } 5
   045 {    45 10  0 } RCL ÷ 0
   046 {       43 44 } g INT
   047 {           2 } 2
   048 {          40 } +
   049 {          36 } ENTER
   050 {          40 } +
   051 {       44  1 } STO 1
   052 {       44  2 } STO 2
   053 {       45  0 } RCL 0
   054 {          16 } CHS
   055 {       44  0 } STO 0
   056 {       43 35 } g CLx
   057 {    42 21  1 } f LBL 1
   058 {       45  2 } RCL 2
   059 {       45  0 } RCL 0
   060 {          14 } y^x
   061 {          30 } −
   062 {          16 } CHS
   063 {    42  5  2 } f DSE 2
   064 {       22  1 } GTO 1
   065 {           2 } 2
   066 {    45 20  0 } RCL × 0
   067 {           1 } 1
   068 {          30 } −
   069 {           2 } 2
   070 {           4 } 4
   071 {          10 } ÷
   072 {       45  1 } RCL 1
   073 {       43 11 } g x²
   074 {          10 } ÷
   075 {           1 } 1
   076 {    45 30  0 } RCL − 0
   077 {           8 } 8
   078 {          10 } ÷
   079 {    45 10  1 } RCL ÷ 1
   080 {          40 } +
   081 {          48 } .
   082 {           5 } 5
   083 {          40 } +
   084 {    45 40  1 } RCL + 1
   085 {       45  0 } RCL 0
   086 {          14 } y^x
   087 {           2 } 2
   088 {          10 } ÷
   089 {          40 } +
   090 {           2 } 2
   091 {       45  0 } RCL 0
   092 {          16 } CHS
   093 {           1 } 1
   094 {          30 } −
   095 {          14 } y^x
   096 {           1 } 1
   097 {          30 } −
   098 {          10 } ÷
   099 {          40 } +

# --------------------------------------------

Hint: replaced the strange "+" and "x" characters generated by the emulator with more common ones.

Hint2: the emulator's extended precision is great, but it seems to use binary arithmetics: enter x=1,05 and the x>y? test that compares x–1 with 0,05 tests true (!) and jumps to LBL 0.

Dieter