Riemann's Zeta Function - another approach (RPL)

[Gerson W. Barbosa]

(06-25-2017 01:16 PM)Dieter Wrote:  

(06-22-2017 01:20 AM)Gerson W. Barbosa Wrote:  For 1 < x < 1.001:

I assume this is 1 < x < 1,01.

(06-22-2017 01:20 AM)Gerson W. Barbosa Wrote:  1.009999999 GSB C -> 100.5779539
1.003456789 GSB C -> 289.8632757

I think we can top this. ;-)

First of all, replace the 13,73328 with 13,7418. Evaluated exactly, this should keep the error within about ±0,7 ULP of a ten-digit result.

But we can do even better:

Zeta(x)  ~  1/u + u/(u + 13,733) + 0,57721568
where u = x–1 and 1 < x ≤ 1,01

Code:

01 LBL C
02 FRAC
03 1/x
04 LastX
05 LastX
06 13,733
07 +
08 /
09 ,57721568
10 +
11 +
12 RTN

1,009999999 => 100,5779533
1,009876543 => 101,8279365
1,003456789 => 289,8632756

According to some quick-and-dirty tests, within the given domain and with exact evaluation (!) the result should have ten significant digits ±0,2 ULP.

If you prefer the result to be either correctly rounded or truncated (error always negative) replace 0,57721568 with ...66:

Zeta (1,001001001) = 999,5772895 48

1,001001001=> 999,5772896 with ...68
1,001001001=> 999,5772895 with ...66

If that's still too much, try
Zeta(x) ~ 1/u + u/(0,9 · u + 13,7335) + 0,577215668

(resp. ...666, see above)

Here the largest error should be about ±0,02 ULP. If evaluated exactly, that is. But this improvement does not always show up on a ten-digit calculator. Actually the limiting factor seems to be the built-in 1/x function with its inherent error of 0,5 ULP.

Thank you for your analysis and improvement. That's what I'll use if I ever try this on the HP-41C.

On the HP 50g I have used 5 Stieltjes constants for 1 < x <= 1.3. Since the size of the numbers doesn't matter much I've kept all constants in the equivalent Horner expression with twelve significant digits.

Code:


%%HP: T(3)A(D)F(.);
DIR
  ZetaX
  \<< RCLF SWAP RAD DUP NOT { .000000000001 + } IFT DUP .5 < { DUP \pi * 2. / SIN DUP + 1. ROT - SWAP OVER GAMMA * \pi DUP + PICK3 ^ / SWAP Zeta * } { Zeta } IFTE SWAP STOF
  \>>
  Zeta
  \<< DUP DUP 1. - ABS .3 > { -1.27 ^ 157. * 1. + 2. / IP DUP + DUP2 0. SWAP 2. / 1.
    FOR k OVER k DUP + DUP 1. - PICK3 ^ INV UNROT SWAP ^ INV SWAP UNROT - + -1.
    STEP 4. ROLLD 3. * 1. + OVER 3. ^ 30. * / PICK3 DUP + 1. + PICK3 SQ 24. * / - PICK3 1. + PICK3 8. * / + .5 + + OVER ^ DUP + INV ROT + 2. ROT ^ DUP 2. - / * } { DUP DUP2 DUP2 -6.66328326918E-6 * 1.3020683574E-4 + * 7.96500246987E-4 - * 3.17028903723E-3 - * 8.10584133725E-2 + * .500000497261 + SWAP 1. - INV + NIP } IFTE
  \>>
END

Zeta: x > 1;

ZetaX: x <> 1.

Real arguments only.

We can use this to check whether 1 + 2 + 3 + 4 + 5 + ... = -1/12    :-)     (YouTube video)

1 +/- ZetaX -> -8.33333333338E-2

1/x -> -11.9999999999

Gerson.