Riemann's Zeta Function - another approach (RPL)

[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