Riemann's Zeta Function - another approach (RPL)

[Dieter]

(07-10-2017 11:49 AM)Dieter Wrote:  It's not perfect until it's perfect. ;-)

Finally I wanted to find out what's possible on a 10-digit calculator, e.g. the HP-41 series. An optimized approximation alone is not sufficient, also some intermediate results have to be protected against roundoff, and even the calculation of 1/(x–1) is not trivial. So the following program tries to carry at least 11 digits if possible, and for 1/(x–1) the integer and fractional parts are evaluated separately.

Code:

 01 LBL "ZETA"       ; Calculates Zeta(x) for 0 ≤ x ≤ 1.106212299
 02 1
 03 -
 04 ENTER
 05 ENTER
 06 ENTER
 07 -8.0516 E-7
 08 *
 09 7.35598 E-6
 10 -
 11 *
 12 9.63289 E-5
 13 +
 14 *
 15 3.4251354 E-4
 16 -
 17 *
 18 4.84521318 E-3
 19 -
 20 *
 21 7.281584515 E-2
 22 +
 23 *
 24 7.721566507 E-2  ; coefficient c0 actually is 0.57721566507 (11 digits)
 25 +                ; the missing 0.5 is added later.
 26 STO 00           ; this yields one more digit for the polynomial.
 27 RDN
 28 1/X
 29 INT
 30 ENTER
 31 RDN
 32 *
 33 CHS
 34 1
 35 +
 36 X<>Y
 37 /
 38 RCL 00
 39 +
 40 X<>Y
 41 .5
 42 +
 43 +
 44 END

The program requires one data register. I am sure with some tricky stack manipulation it can be done without it, but the way it is the program can be used on most classic HPs (after you have prestored the constants somewhere).

The implemented approximation is good for an error within ±6 units in the 12th significant digit, so with sufficient precision the rounded 10-digit result should be within 0,56 ULP.

After some first tests it looks like the returned 10-digit result indeed is within about ±0,6 ULP, i.e. comparable to hardware accuracy. Here and there the result may be 1 ULP off, e.g. at Zeta(0,01). The true 12-digit value is -0.509290714040 while the program calculates an 11-digit intermediate result of ....71405, so the final result is rounded up to ...41 instead of down to ...40. #-) All in all the returned values look very good, I think. If you find any larger errors please report here.

Dieter