Riemann's Zeta Function - another approach (RPL)

[Gerson W. Barbosa]

Three lines shorter:

Code:


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

   000 {             } 
   001 {    42 21 11 } f LBL A
   002 {       44  3 } STO 3
   003 {           1 } 1
   004 {          30 } -
   005 {       43 16 } g ABS
   006 {          48 } .
   007 {           0 } 0
   008 {           4 } 4
   009 {          34 } x↔y
   010 {       43 10 } g x≤y
   011 {       22  1 } GTO 1
   012 {       45  3 } RCL 3
   013 {          36 } ENTER
   014 {          36 } ENTER
   015 {          48 } .
   016 {           8 } 8
   017 {           2 } 2
   018 {          16 } CHS
   019 {          14 } y^x
   020 {           2 } 2
   021 {           5 } 5
   022 {          20 } ×
   023 {       43 44 } g INT
   024 {          36 } ENTER
   025 {          40 } +
   026 {       44 25 } STO I
   027 {       44  2 } STO 2
   028 {          34 } x↔y
   029 {          16 } CHS
   030 {       44  1 } STO 1
   031 {           0 } 0
   032 {       44  0 } STO 0
   033 {    42 21  0 } f LBL 0
   034 {       45 25 } RCL I
   035 {       45  1 } RCL 1
   036 {          14 } y^x
   037 {    44 30  0 } STO - 0
   038 {           1 } 1
   039 {    44 30 25 } STO - I
   040 {       45 25 } RCL I
   041 {       45  1 } RCL 1
   042 {          14 } y^x
   043 {    44 40  0 } STO + 0
   044 {    42  5 25 } f DSE I
   045 {       22  0 } GTO 0
   046 {       45  1 } RCL 1
   047 {          36 } ENTER
   048 {          40 } +
   049 {           1 } 1
   050 {          30 } -
   051 {       45  2 } RCL 2
   052 {       43 11 } g x²
   053 {           2 } 2
   054 {           4 } 4
   055 {          20 } ×
   056 {          10 } ÷
   057 {           1 } 1
   058 {    45 30  1 } RCL - 1
   059 {           8 } 8
   060 {    45 20  2 } RCL × 2
   061 {          10 } ÷
   062 {          40 } +
   063 {          48 } .
   064 {           5 } 5
   065 {          40 } +
   066 {    45 40  2 } RCL + 2
   067 {       45  1 } RCL 1
   068 {          14 } y^x
   069 {           2 } 2
   070 {          10 } ÷
   071 {    45 40  0 } RCL + 0
   072 {           2 } 2
   073 {       45  1 } RCL 1
   074 {          16 } CHS
   075 {          14 } y^x
   076 {          36 } ENTER
   077 {          36 } ENTER
   078 {           2 } 2
   079 {          30 } -
   080 {          10 } ÷
   081 {          20 } ×
   082 {       43 32 } g RTN
   083 {    42 21  1 } f LBL 1
   084 {       45  3 } RCL 3
   085 {           1 } 1
   086 {          30 } -
   087 {          15 } 1/x
   088 {       43 36 } g LSTx
   089 {          48 } .
   090 {           9 } 9
   091 {       43 36 } g LSTx
   092 {          20 } ×
   093 {           1 } 1
   094 {           3 } 3
   095 {          48 } .
   096 {           7 } 7
   097 {           3 } 3
   098 {           3 } 3
   099 {           4 } 4
   100 {           4 } 4
   101 {          40 } +
   102 {          10 } ÷
   103 {          48 } .
   104 {           5 } 5
   105 {           7 } 7
   106 {           7 } 7
   107 {           2 } 2
   108 {           1 } 1
   109 {           5 } 5
   110 {           6 } 6
   111 {           7 } 7
   112 {          40 } +
   113 {          40 } +
   114 {       43 32 } g RTN
   115 {    42 21 12 } f LBL B
   116 {          48 } .
   117 {           5 } 5
   118 {          34 } x↔y
   119 {    43 30  0 } g TEST x≠0
   120 {       22  2 } GTO 2
   121 {          34 } x↔y
   122 {          16 } CHS
   123 {       43 32 } g RTN
   124 {    42 21  2 } f LBL 2
   125 {       43 10 } g x≤y
   126 {       22  3 } GTO 3
   127 {       32 11 } GSB A
   128 {       43 32 } g RTN
   129 {    42 21  3 } f LBL 3
   130 {           1 } 1
   131 {          34 } x↔y
   132 {          30 } -
   133 {       44  4 } STO 4
   134 {       32 11 } GSB A
   135 {       43 26 } g π
   136 {          36 } ENTER
   137 {          40 } +
   138 {       45  4 } RCL 4
   139 {          14 } y^x
   140 {          10 } ÷
   141 {           1 } 1
   142 {    45 30  4 } RCL - 4
   143 {           1 } 1
   144 {       43 23 } g SIN-¹
   145 {          20 } ×
   146 {          23 } SIN
   147 {          20 } ×
   148 {           1 } 1
   149 {          16 } CHS
   150 {    45 40  4 } RCL + 4
   151 {       42  0 } f x!
   152 {          20 } ×
   153 {          36 } ENTER
   154 {          40 } +
   155 {       43 32 } g RTN

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

\[\zeta(x)=\frac{2\cdot (-x)!\cdot \sin \left ( x\cdot \sin^{-1} 1 \right )\zeta (1-x)}{(2 \pi )^{1-x}}\]

[Dieter]

(06-28-2017 07:45 PM)Gerson W. Barbosa Wrote:  Three lines shorter:

Here is my first attempt at a HP41 version for x > 1. Maybe this or that can be reused in your 15C program, e.g. the use of less registers (the 15C would require R0...R2) or the DSE within the loop.

Code:

01  LBL"ZETA"
02  STO 00
03  1
04  -
05  LN     // generate error if x<=1
06  ,01
07  LastX
08  x>y?
09  GTO 00
10  1/x
11  LastX
12  LastX
13  ,9
14  *
15  13,7335
16  +
17  /
18  ,57721567
19  +
20  +
21  GTO 02
22  LBL 00
23  28
24  RCL 00
25  /
26  INT
27  1
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  DSE Y
41  RCL Y
42  RCL 00
43  y^x
44  +
45  DSE Y
46  GTO 01
47  RCL 00
48  ST+ X
49  1
50  -
51  24
52  /
53  RCL 01
54  x^2
55  /
56  1
57  RCL 00
58  -
59  8
60  /
61  RCL 01
62  /
63  +
64  ,5
65  +
66  RCL 01
67  +
68  RCL 00
69  y^x
70  2
71  /
72  +
73  1
74  RCL 00
75  +
76  2
77  LN
78  *
79  e^x-1
80  CHS
81  /
82  LBL 02
83  END

Edit: corrected an error in line 47

Your method for calculating the number of required terms of the sum is fine, actually I got the same result. However the program uses a shorter formula that yields very similar results. The final calculation with this power-of-2-fraction can cause significant errors for x close to 1, so I chose a different approach with the 41's e^x–1 function which should reduce the error here. The simplified Zeta evaluation is used for 1 < x ≤ 1,01. At x=1,01 both methods agree and return the correct result 100,5779433.

The examples in your first post (except the two cases that cannot be handled by the program above) are all evaluated to ten correct digits. But of course the program is not always that accurate. ;-)

Dieter

[Gerson W. Barbosa]

(06-29-2017 07:47 PM)Dieter Wrote:  Here is my first attempt at a HP41 version for x > 1. Maybe this or that can be reused in your 15C program, e.g. the use of less registers (the 15C would require R0...R2) or the DSE within the loop.

Code:

01  LBL"ZETA"
02  STO 00
03  1
04  -
05  LN     // generate error if x<=1
06  ,01
07  LastX
08  x>y?
09  GTO 00
10  1/x
11  LastX
12  LastX
13  ,9
14  *
15  13,7335
16  +
17  /
18  ,57721567
19  +
20  +
21  GTO 02
22  LBL 00
23  28
24  RCL 00
25  /
26  INT
27  1
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  DSE Y
41  RCL Y
42  RCL 00
43  y^x
44  +
45  DSE Y
46  GTO 01
47  RCL 01
48  ST+ X
49  1
50  -
51  24
52  /
53  RCL 01
54  x^2
55  /
56  1
57  RCL 00
58  -
59  8
60  /
61  RCL 01
62  /
63  +
64  ,5
65  +
66  RCL 01
67  +
68  RCL 00
69  y^x
70  2
71  /
72  +
73  1
74  RCL 00
75  +
76  2
77  LN
78  *
79  e^x-1
80  CHS
81  /
82  LBL 02
83  END

Your method for calculating the number of required terms of the sum is fine, actually I got the same result. However the program uses a shorter formula that yields very similar results. The final calculation with this power-of-2-fraction can cause significant errors for x close to 1, so I chose a different approach with the 41's e^x–1 function which should reduce the error here. The simplified Zeta evaluation is used for 1 < x ≤ 1,01. At x=1,01 both methods agree and return the correct result 100,5779433.

The examples in your first post (except the two cases that cannot be handled by the program above) are all evaluated to ten correct digits. But of course the program is not always that accurate. ;-)

Both your method and mine causes the series to be calculated to more terms than are really needed. For instance, Zeta(2) requires only 22 terms, but we evaluate 28 and 30 terms, respectively. A closer fit would require a longer formula, so yours presents a good compromise between size and exactness.

1.2 40
2 22
3 18
4 16
5 12
6 10
7 8
8 8
9 6
10 6
...
12 6
13 4
...
20 4
21 2
...
40 2

For the 15C I gave two programs, A and B. If there's no need for the analytical continuation then B should be ignored. But A should be fast enough for x >= 1/2 if it is supposed to be called by B. Your HP-41 programs computes Zeta(1/2) in about one minute, which is pretty good.

There is a typo in your line 47: RCL 00, not 01.

Thank you for providing a 41 version. I will try it on the 42S.

Gerson.

[Dieter]

(06-30-2017 02:20 AM)Gerson W. Barbosa Wrote:  Both your method and mine causes the series to be calculated to more terms than are really needed. For instance, Zeta(2) requires only 22 terms, but we evaluate 28 and 30 terms, respectively. A closer fit would require a longer formula, so yours presents a good compromise between size and exactness.

Yes, especially between 1 and 2,5 the calculated number of terms is higher than required. BTW, how did you determine the exact numbers here?

(06-30-2017 02:20 AM)Gerson W. Barbosa Wrote:  There is a typo in your line 47: RCL 00, not 01.

Yes, thank you. This has been corrected now.

Dieter

[Gerson W. Barbosa]

(06-30-2017 06:42 AM)Dieter Wrote:  

(06-30-2017 02:20 AM)Gerson W. Barbosa Wrote:  Both your method and mine causes the series to be calculated to more terms than are really needed. For instance, Zeta(2) requires only 22 terms, but we evaluate 28 and 30 terms, respectively. A closer fit would require a longer formula, so yours presents a good compromise between size and exactness.

Yes, especially between 1 and 2,5 the calculated number of terms is higher than required. BTW, how did you determine the exact numbers here?

I tested one by one with ever higher n until I got all 10 significant digits right. I did this on the 15C LE which is 150x faster than the original one. I could have used the HP-15C iOS app as it is even faster:

0.0001 GSB A -> -0.5000919(117) (8 s, 95272 terms)

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.

Thanks again for your great contributions, especially for your striving to always get as many significant digits as possible.

Gerson.

[Paul Dale]

How does this compare to Jean-Marc Baillard's implementation of Borwein's second algorithm?

The 34S uses this algorithm. Originally in C but later in XROM. Borwein's paper includes an error term which means that for real arguments, the number of terms for a specified precision is constant & can be determined in advance. This isn't true for complex numbers, where the number of terms depends on the magnitude of the complex part.


Pauli

[Gerson W. Barbosa]

(06-30-2017 12:08 PM)Paul Dale Wrote:  How does this compare to Jean-Marc Baillard's implementation of Borwein's second algorithm?

Quoting from your first link:

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

   3  XEQ "ZETA" ->  1.202056903    (15 s) [HP-41CV]
-7.49 GSB    B   ->  0.003312040168 (26 s) [HP-15C] (probably 11 seconds on the HP-41CV)
 1.1  XEQ "ZETA" -> 10.58444846     (34 s) [HP-41CV]


This relies on an empirical correction expression I've found, though:

1/(2*((n + 1/2 + (s + 1)/(8*n) - (2*s + 1)/(24(?)*n^2) + ... )^x))

But I am still not sure wheter the last term is correct or if this is correct at all...

Gerson.

PS: Perhaps Borwein's algorithm is overkill for the HP-41. If more digits are to be calculated, as on the wp34s, then is should be faster, even if more terms of the correction expression were available.

[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

[Paul Dale]

(06-30-2017 04:50 PM)Dieter Wrote:  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)

How about \( e^x - 1 = tanh\left(\frac{x}{2}\right) \left( e^x + 1 \right) \)?

Pauli

[Dieter]

(06-30-2017 11:23 PM)Paul Dale Wrote:  How about \( e^x - 1 = tanh\left(\frac{x}{2}\right) \left( e^x + 1 \right) \)?

Ah, this reminds me of an earlier thread on e^x–1.

In this special case the function may be coded this way, cf. Gerson's post #19 in the mentioned thread:

Code:

 .. ...         ...
 90 2           1
 91 RCL 0       STO+ 0
 92 1           2
 93 +           RCL 0
 94 y^x         y^x
 95 LastX       +
 96 2           /
 97 LN          2
 98 *           LN
 99 2           2
100 /           /
101 TANH        RCL* 0
102 *           TANH
103 Last        /
104 +           CHS
105 /
106 CHS

Do not use ln √2 instead of ½ ln 2 here as on the 10 digit calculators this causes an error of several ULPs.

Dieter

[Gerson W. Barbosa]

(06-30-2017 04:50 PM)Dieter Wrote:  

(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?

Oops! I forgot the HP-41 lacked Gamma.

(06-30-2017 04:50 PM)Dieter Wrote:  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.

I chose to evaluate two terms per loop because it's faster. For example, your older program takes 65 seconds to compute Zeta(1/2), if your restriction is by-passed, while your new one takes 73 seconds. Instead of the first DSE, I tried SIGN ST+Y X<>L, but this would take 70 seconds, so I kept DSE. It's not a true DEC (decrement) instruction as it includes an implicit time consuming IF, but it's faster than what I'd been using in the 15C program.

(06-30-2017 04:50 PM)Dieter Wrote:  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

I have used it in my new version. Not all your suggestions have been implemented yet.

Code:


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

   000 {             } 
   001 {    42 21 11 } f LBL A
   002 {          48 } .
   003 {           5 } 5
   004 {          34 } x↔y
   005 {    43 30  0 } g TEST x≠0
   006 {       22  2 } GTO 2
   007 {          34 } x↔y
   008 {          16 } CHS
   009 {       43 32 } g RTN
   010 {    42 21  2 } f LBL 2
   011 {       43 10 } g x≤y
   012 {       22  3 } GTO 3
   013 {       32  4 } GSB 4
   014 {       43 32 } g RTN
   015 {    42 21  3 } f LBL 3
   016 {           1 } 1
   017 {          34 } x↔y
   018 {          30 } -
   019 {       44  4 } STO 4
   020 {       32  4 } GSB 4
   021 {       43 26 } g π
   022 {          36 } ENTER
   023 {          40 } +
   024 {       45  4 } RCL 4
   025 {          14 } y^x
   026 {          10 } ÷
   027 {           1 } 1
   028 {    45 30  4 } RCL - 4
   029 {           1 } 1
   030 {       43 23 } g SIN-¹
   031 {          20 } ×
   032 {          23 } SIN
   033 {          20 } ×
   034 {           1 } 1
   035 {          16 } CHS
   036 {    45 40  4 } RCL + 4
   037 {       42  0 } f x!
   038 {          20 } ×
   039 {          36 } ENTER
   040 {          40 } +
   041 {       43 32 } g RTN
   042 {    42 21  4 } f LBL 4
   043 {       44  3 } STO 3
   044 {           1 } 1
   045 {          30 } -
   046 {       43 16 } g ABS
   047 {          48 } .
   048 {           0 } 0
   049 {           5 } 5
   050 {          34 } x↔y
   051 {       43 10 } g x≤y
   052 {       22  1 } GTO 1
   053 {       45  3 } RCL 3
   054 {          36 } ENTER
   055 {          36 } ENTER
   056 {          48 } .
   057 {           8 } 8
   058 {           2 } 2
   059 {          16 } CHS
   060 {          14 } y^x
   061 {           2 } 2
   062 {           5 } 5
   063 {          20 } ×
   064 {       43 44 } g INT
   065 {          36 } ENTER
   066 {          40 } +
   067 {       44 25 } STO I
   068 {       44  2 } STO 2
   069 {          34 } x↔y
   070 {          16 } CHS
   071 {       44  1 } STO 1
   072 {           0 } 0
   073 {       44  0 } STO 0
   074 {    42 21  0 } f LBL 0
   075 {       45 25 } RCL I
   076 {       45  1 } RCL 1
   077 {          14 } y^x
   078 {    44 30  0 } STO - 0
   079 {    42  5 25 } f DSE I
   080 {       45 25 } RCL I
   081 {       45  1 } RCL 1
   082 {          14 } y^x
   083 {    44 40  0 } STO + 0
   084 {    42  5 25 } f DSE I
   085 {       22  0 } GTO 0
   086 {       45  1 } RCL 1
   087 {          36 } ENTER
   088 {          40 } +
   089 {           1 } 1
   090 {          30 } -
   091 {       45  2 } RCL 2
   092 {       43 11 } g x²
   093 {           2 } 2
   094 {           4 } 4
   095 {          20 } ×
   096 {          10 } ÷
   097 {           1 } 1
   098 {    45 30  1 } RCL - 1
   099 {           8 } 8
   100 {    45 20  2 } RCL × 2
   101 {          10 } ÷
   102 {          40 } +
   103 {          48 } .
   104 {           5 } 5
   105 {          40 } +
   106 {    45 40  2 } RCL + 2
   107 {       45  1 } RCL 1
   108 {          14 } y^x
   109 {           2 } 2
   110 {          10 } ÷
   111 {    45 40  0 } RCL + 0
   112 {           2 } 2
   113 {       45  1 } RCL 1
   114 {          16 } CHS
   115 {           1 } 1
   116 {          30 } -
   117 {          14 } y^x
   118 {           1 } 1
   119 {          30 } -
   120 {          10 } ÷
   121 {          40 } +
   122 {       43 32 } g RTN
   123 {    42 21  1 } f LBL 1
   124 {       45  3 } RCL 3
   125 {           1 } 1
   126 {          30 } -
   127 {          15 } 1/x
   128 {       43 36 } g LSTx
   129 {          48 } .
   130 {           9 } 9
   131 {           2 } 2
   132 {           4 } 4
   133 {           6 } 6
   134 {       43 36 } g LSTx
   135 {          20 } ×
   136 {           1 } 1
   137 {           3 } 3
   138 {          48 } .
   139 {           7 } 7
   140 {           3 } 3
   141 {           3 } 3
   142 {          40 } +
   143 {          10 } ÷
   144 {          48 } .
   145 {           5 } 5
   146 {           7 } 7
   147 {           7 } 7
   148 {           2 } 2
   149 {           1 } 1
   150 {           5 } 5
   151 {           6 } 6
   152 {           5 } 5
   153 {           4 } 4
   154 {          40 } +
   155 {          40 } +
   156 {       43 32 } g RTN

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

 -1.5  GSB A   ->   -0.02548520190 [ 89] (  47 s)
 -1    GSB A   ->   -0.08333333332 [  3] (  52 s)
 0     GSB A   ->   -0.500000000         (   1 s)
 0.5   GSB A   ->   -1.460354507   [  9] ( 144 s)
 0.02  GSB A   ->   -0.5187882122  [ 14] (  11 s)
 0.94  GSB A   ->  -16.09383730    [  2] (  86 s)
 0.959 GSB A   ->  -23.81602202    [181] (   5 s) *
 0.961 GSB A   ->  -25.06665734    [ 14] (   5 s)
 1.02  GSB A   ->   50.57867004          (   5 s)
 1.04  GSB A   ->   25.58012052          (   5 s)
 1.041 GSB A   ->   24.97043684    [  5] (   5 s)
 1.06  GSB A   ->   17.24823369    [  7] (  83 s)
 2     GSB A   ->    1.644934067         (  48 s)
 3     GSB A   ->    1.202056903         (  36 s)
 4     GSB A   ->    1.082323234         (  31 s)
10     GSB A   ->    1.000994575         (  17 s)

The angle mode is irrelevant.

* Somewhat off! Will see that later. Vacation time mode on :-)

[Dieter]

(07-01-2017 02:57 PM)Gerson W. Barbosa Wrote:  Oops! I forgot the HP-41 lacked Gamma.

That's an essential point here.

(07-01-2017 02:57 PM)Gerson W. Barbosa Wrote:  I chose to evaluate two terms per loop because it's faster. For example, your older program takes 65 seconds to compute Zeta(1/2), if your restriction is by-passed, while your new one takes 73 seconds.

Yes, evaluating multiple terms per loop is faster, especially since the GTO and the required label search take some time. On the 41 this is less of an issue because of the "compiled GTO" feature. But anyway it still makes a small difference.

(07-01-2017 02:57 PM)Gerson W. Barbosa Wrote:  Instead of the first DSE, I tried SIGN ST+Y X<>L, but this would take 70 seconds, so I kept DSE. It's not a true DEC (decrement) instruction as it includes an implicit time consuming IF, but it's faster than what I'd been using in the 15C program.

DSE usually is very fast, not only on the 41. Sure, the implicit test is not required, but all is done in a single command and without stack interference, so this is the preferred choice. For (integer) 2x–1 I often use ST+X DSE X.

Code:

   010 {    42 21  2 } f LBL 2
   011 {       43 10 } g x≤y
   012 {       22  3 } GTO 3
   013 {       32  4 } GSB 4
   014 {       43 32 } g RTN
   015 {    42 21  3 } f LBL 3

This can be recplaced by

LBL 2
X>Y?
GTO 4

Since LBL 3 is not used elsewhere it can be omitted.

Code:

   042 {    42 21  4 } f LBL 4
   043 {       44  3 } STO 3
   044 {           1 } 1
   045 {          30 } -
   046 {       43 16 } g ABS
   047 {          48 } .
   048 {           0 } 0
   049 {           5 } 5
   050 {          34 } x↔y
   051 {       43 10 } g x≤y

This tests whether x is in 0,95...1,05. My approximation is valid for 1...1,05. It hasn't been designed for x<1. Here the error is substantially larger, that's why your results for x=0,959 and 0,961 are a bit off in the last digits. At this point the error is about 20 ULP, and that's exactly what your results show.

Code:

   053 {       45  3 } RCL 3
   054 {          36 } ENTER
   055 {          36 } ENTER
   056 {          48 } .
   057 {           8 } 8
   058 {           2 } 2
   059 {          16 } CHS
   060 {          14 } y^x
   061 {           2 } 2
   062 {           5 } 5
   063 {          20 } ×
   064 {       43 44 } g INT
   065 {          36 } ENTER
   066 {          40 } +
   067 {       44 25 } STO I
   068 {       44  2 } STO 2
   069 {          34 } x↔y
   070 {          16 } CHS
   071 {       44  1 } STO 1

First, the calculated number of terms has no added constant, so for x>50,7 it will evaluate to zero.
Then, why do you first pust R3 onto the stack and then pop it back later?

2
5
RCL 3
,
8
2
y^x
/
INT
1
+
ENTER
+
STO I
STO 2
RCL 3
CHS
STO 1

Also another register for the summed terms is not required, you can simply do it on the stack. ;-)

(07-01-2017 02:57 PM)Gerson W. Barbosa Wrote:  * Somewhat off! Will see that later.

cf. my comment above.

(07-01-2017 02:57 PM)Gerson W. Barbosa Wrote:  Vacation time mode on :-)

In this case: have a good time. :-)

Dieter

[Gerson W. Barbosa]

Thank you very much, Dieter, for kind of doing my homework:-)

l'm looking foward for other implementations which might demonstrate whether this approach is worthwhile or not. I wasn't aware of J. M. Baillard's new programs, else probably I would't even have tried all this.

Gerson

[Dieter]

(07-02-2017 12:52 PM)Gerson W. Barbosa Wrote:  Thank you very much, Dieter, for kind of doing my homework:-)

I also did mine, and so here is a new HP41 program.

First, there is a new approximation for 0,97≤x≤1,03 with an error of approx. ±0,5 units in the 10th significant digit.

Then I realized that for 0,3454<x<0,97 the program may use a constant number of iterations and the results show a relatively constant error (only a few ULP) that can be compensated by a heuristc formula. Here I chose 54 terms, and the result is corrected by an amount of approx. 5 E–9/x^¼. This keeps the error within about ±0,6 ULP.

The lower limit for x is the point where Zeta(x) is exactly –1. Beyond this the accuracy will substantially degrade, so x=0,3453726573 is the limit here. Below this the program throws a DATA ERROR.

So here is the latest version:

Code:

 01 LBL "ZETA"
 02 STO 00
 03 .3453726573
 04 -
 05 SQRT
 06 .03
 07 RCL 00
 08 1
 09 -
 10 ABS
 11 X>Y?
 12 GTO 00
 13 LASTX
 14 1/X
 15 LASTX
 16 LASTX
 17 .9135
 18 *
 19 13.73336
 20 +
 21 /
 22 .577215664
 23 +
 24 +
 25 GTO 02
 26 LBL 00
 27 26
 28 RCL 00
 29 /
 30 2
 31 +
 32 INT
 33 ST+ X
 34 54
 35 X>Y?
 36 X<>Y
 37 STO 01
 38 RCL 00
 39 CHS
 40 STO 00
 41 CLX
 42 LBL 01
 43 RCL Y
 44 RCL 00
 45 Y^X
 46 -
 47 CHS
 48 DSE Y
 49 GTO 01
 50 RCL 00
 51 ST+ X
 52 1
 53 -
 54 24
 55 /
 56 RCL 01
 57 X^2
 58 /
 59 1
 60 RCL 00
 61 -
 62 8
 63 /
 64 RCL 01
 65 /
 66 +
 67 .5
 68 +
 69 RCL 01
 70 +
 71 RCL 00
 72 Y^X
 73 2
 74 /
 75 +
 76 RCL 00
 77 1
 78 +
 79 2
 80 LN
 81 *
 82 E^X-1
 83 CHS
 84 /
 85 RCL 00
 86 1
 87 +
 88 SIGN
 89 X<0?
 90 ST- X
 91 5.5 E-9
 92 *
 93 RCL 00
 94 ABS
 95 SQRT
 96 SQRT
 97 /
 98 -
 99 LBL 02
100 END

Examples:

Code:

0,35  XEQ"ZETA"  => -1,010511224   exact   37 s
0,5   XEQ"ZETA"  => -1,460354509   exact   37 s
0,75  XEQ"ZETA"  => -3,441285386  -1 ULP   37 s
0,9   XEQ"ZETA"  => -9,430114018  -1 ULP   37 s
0,97  XEQ"ZETA"  => -32,75830650   exact    2 s
0,999 XEQ"ZETA"  => -999,4228572   exact    2 s
1,001 XEQ"ZETA"  =>  1000,577288   exact    2 s
1,03  XEQ"ZETA"  =>  33,91272911  +1 ULP    2 s
1,05  XEQ"ZETA"  =>  20,58084429  -1 ULP   36 s
1,27  XEQ"ZETA"  =>  4,300220201   exact   31 s
2     XEQ"ZETA"  =>  1,644934067   exact   22 s
3     XEQ"ZETA"  =>  1,202056903   exact   16 s
5     XEQ"ZETA"  =>  1,036927755   exact   12 s
19,99 XEQ"ZETA"  =>  1,000000961   exact    8 s
30    XEQ"ZETA"  =>  1,000000001   exact    6 s

Finally Zeta(0,3453726573) returns exactly –1.

The given execution times are due to V41 and a speed setting that quite exactly matches the real thing.

Dieter

[Gerson W. Barbosa]

(07-02-2017 02:54 PM)Dieter Wrote:  First, there is a new approximation for 0,97≤x≤1,03 with an error of approx. ±0,5 units in the 10th significant digit.

Then I realized that for 0,3454<x<0,97 the program may use a constant number of iterations and the results show a relatively constant error (only a few ULP) that can be compensated by a heuristc formula. Here I chose 54 terms, and the result is corrected by an amount of approx. 5 E–9/x^¼. This keeps the error within about ±0,6 ULP.

The lower limit for x is the point where Zeta(x) is exactly –1. Beyond this the accuracy will substantially degrade, so x=0,3453726573 is the limit here. Below this the program throws a DATA ERROR.

Great! It needs to be accurate only for x > 0.5. Then the functional equation can handle the rest if Gamma is available.

Your code can be nicely pasted into Free42, except line 91 (because the space between the last digit and E). LBL 2 is not necessary, unless there's a reason to keep it:

Code:


00 { 177-Byte Prgm }
01▸LBL "ZETA"
02 STO 00
03 0,3453726573
04 -
05 SQRT
06 0,03
07 RCL 00
08 1
09 -
10 ABS
11 X>Y?
12 GTO 00
13 LASTX
14 1/X
15 LASTX
16 LASTX
17 0,9135
18 ×
19 13,73336
20 +
21 ÷
22 0,577215664
23 +
24 +
25 RTN
26▸LBL 00
27 26
28 RCL 00
29 ÷
30 2
31 +
32 IP
33 STO+ ST X
34 54
35 X>Y?
36 X<>Y
37 STO 01
38 RCL 00
39 +/-
40 STO 00
41 CLX
42▸LBL 01
43 RCL ST Y
44 RCL 00
45 Y↑X
46 -
47 +/-
48 DSE ST Y
49 GTO 01
50 RCL 00
51 STO+ ST X
52 1
53 -
54 24
55 ÷
56 RCL 01
57 X↑2
58 ÷
59 1
60 RCL 00
61 -
62 8
63 ÷
64 RCL 01
65 ÷
66 +
67 0,5
68 +
69 RCL 01
70 +
71 RCL 00
72 Y↑X
73 2
74 ÷
75 +
76 RCL 00
77 1
78 +
79 2
80 LN
81 ×
82 E↑X-1
83 +/-
84 ÷
85 RCL 00
86 1
87 +
88 SIGN
89 X<0?
90 STO- ST X
91 5,5ᴇ-9
92 ×
93 RCL 00
94 ABS
95 SQRT
96 SQRT
97 ÷
98 -
99 END

Gerson.

[Dieter]

(07-02-2017 11:53 PM)Gerson W. Barbosa Wrote:  Great! It needs to be accurate only for x > 0.5. Then the functional equation can handle the rest if Gamma is available.

I have tried also this on the '41, using the FACT function so that at least integer arguments can be handled, and it worked fine. Now, if there only was a Gamma function...

BTW, does anyone know if there was a ROM that included Gamma in Mcode (!) ? Except the ones by Ángel, that is. ;-) The Advantage ROM was said to update the '41 to features that else were only available on the 34C and 15C (Solve, Integrate, Matrix operations), but AFAIK this did not include Gamma.

(07-02-2017 11:53 PM)Gerson W. Barbosa Wrote:  Your code can be nicely pasted into Free42, except line 91 (because the space between the last digit and E).

The listing was generated by hp41uc. But remember that all the math is designed for 10-digit precision. On the 42s you may be better off with a constant 5 E–9. Try it and see what you get. Maybe you want to round the results to 10 digits by SCI 9 RND ALL. ;-) On the 42s you can also apply RCL-arithmetics which will save a few steps.

(07-02-2017 11:53 PM)Gerson W. Barbosa Wrote:  LBL 2 is not necessary, unless there's a reason to keep it:

The reason is a common exit point at the end of the program. This way it can be restarted with a simple R/S. ;-)

Dieter

[Gerson W. Barbosa]

(07-03-2017 07:07 AM)Dieter Wrote:  I have tried also this on the '41, using the FACT function so that at least integer arguments can be handled, and it worked fine. Now, if there only was a Gamma function...

BTW, does anyone know if there was a ROM that included Gamma in Mcode (!) ? Except the ones by Ángel, that is. ;-) The Advantage ROM was said to update the '41 to features that else were only available on the 34C and 15C (Solve, Integrate, Matrix operations), but AFAIK this did not include Gamma.

...

The listing was generated by hp41uc. But remember that all the math is designed for 10-digit precision. On the 42s you may be better off with a constant 5 E–9. Try it and see what you get. Maybe you want to round the results to 10 digits by SCI 9 RND ALL. ;-) On the 42s you can also apply RCL-arithmetics which will save a few steps.

...

(07-02-2017 11:53 PM)Gerson W. Barbosa Wrote:  LBL 2 is not necessary, unless there's a reason to keep it:

The reason is a common exit point at the end of the program. This way it can be restarted with a simple R/S. ;-)

1) Unfortunately Gamma didn't make into the Advantage ROM. That's the only ROM module I have. Really disappointing.
I guess you haven't considered RPN implementations because they are rather slow.

2) Yes, I am aware of that. I used Free42 only to check your new version on what I had at hand.

3) That makes quite sense! Anyway, on the 42 I usually assign the program I am most interested at the moment to the Sigma+ key to avoid the lengthy XEQ ALPHA XXXX ALPHA sequence.

Gerson.

[Dieter]

(07-02-2017 02:54 PM)Dieter Wrote:  First, there is a new approximation for 0,97≤x≤1,03 with an error of approx. ±0,5 units in the 10th significant digit.

Update: I remembered the Zeta evaluation using a polynomial with the Stieltjes' constants, and so the idea was to develop a dedicated polynomial approximation. Here I chose a simple quadratic equation so that only three constants are required. This is what I came up with:

Zeta(x) ~ 1/u + 0,577215665 + 0,07281553 · u – 0,0048452 · u²
where u = x–1  and  0,965 ≤ x ≤ 1,035

The error within the given domain is less than 4 E–9, i.e. here less than 0,4 ULP (if evaluated with sufficient precision, that is).

One more term (i.e. a third order polynomial) improves the results substantially, you can easily handle 0,9 ≤ x ≤ 1,1 with a possible accuracy better than 0,2 units in the 10th significant digit:

Zeta(x) ~ 1/u + 0,5772156632 + 0,072815845 · u – 0,00484404 · u² – 0,0003424 · u³
where u = x–1 and 0,9 ≤ x ≤ 1,1

One more update:

Think big. ;-) I tried a fourth-order polynomial and the results are quite promising. You can go all the way from 1,1 down to 0,5 (!) with an error less than half a unit in the 10th significant digit. Maybe you want to try this:

Zeta(x) ~ 1/u + 0,5772156625 + 0,072815839 · u – 0,004844823 · u² – 0,000339474 · u³ + 0,000104308 · u⁴
where u = x–1 and 0,5 ≤ x ≤ 1,1

Dieter

[Dieter]

(07-03-2017 07:30 PM)Dieter Wrote:  I tried a fourth-order polynomial and the results are quite promising. You can go all the way from 1,1 down to 0,5 (!) with an error less than half a unit in the 10th significant digit.

If you got a calculator with, say, 12 digit precision. ;-)

Anyway, here is a final addition. It is a fifth order approximation for 0≤x≤1,05 which has been tested in Excel with every intermediate result rounded to 10 significant digits. And indeed the results seem to match those on a real 15C or 41C. If the implementation is carefully coded, cf. step 23 ff. in the listing below, the results should be within 2 units in the 10th place. Which I think is as good as it gets with ten-digit working precision.

The formula:

Zeta(x) = 1/u + 0,577215668 + 0,0728159383 · u – 0,0048444781 · u² – 0,0003401389 · u³ + 0,0001000277 · u⁴ – 4,58184E-6 · u⁵
where u = x – 1  and  0 ≤ x ≤ 1,05.

Edit: please note the slightly improved coefficients in post #42.

After storing the constants in R0 (=0,5772...) to R5 (=–4,58184E-6) – be sure to observe the correct signs – the following program should yield results within 2 ULP. If you find cases with larger errors, please report here.

Code:

01 LBL A
02 1
03 -
04 ENTER
05 ENTER
06 ENTER
07 RCL 5
08 *
09 RCL 4
10 +
11 *
12 RCL 3
13 +
14 *
15 RCL 2
16 +
17 *
18 RCL 1
19 +
20 *
21 RCL 0
22 +
23 1
24 R^
25 +
26 LstX
27 /
28 +
29 1
32 -
33 RTN

Examples:

1,05 => 20,58084431  (+1 ULP)
0,5  => -1,460354508 (–1 ULP)
 0   => -0,5000000000 (exact)

Dieter

[Gerson W. Barbosa]

(07-07-2017 07:40 PM)Dieter Wrote:  Anyway, here is a final addition. It is a fifth order approximation for 0≤x≤1,05 which has been tested in Excel with every intermediate result rounded to 10 significant digits. And indeed the results seem to match those on a real 15C or 41C. If the implementation is carefully coded, cf. step 23 ff. in the listing below, the results should be within 2 units in the 10th place. Which I think is as good as it gets with ten-digit working precision.

The formula:

Zeta(x) = 1/u + 0,577215668 + 0,0728159383 · u – 0,0048444781 · u² – 0,0003401389 · u³ + 0,0001000277 · u⁴ – 4,58184E-6 · u⁵
where u = x – 1  and  0 ≤ x ≤ 1,05.

Very nice!

http://m.wolframalpha.com/input/?i=plot+...05&x=2&y=6

Do you have something as good as that for 12-digits calculators? Thanks!

Gerson.